scaling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "899f4141",
   "metadata": {},
   "source": [
    "# Scaling from mendes santos paper\n",
    "1. **import data** \n",
    "1. **define scaling functions**\n",
    "1. **make parameter mesh**\n",
    "1. **plot 3d scatter**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 379,
   "id": "a87c160a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from fssa import quality\n",
    "from tqdm import tqdm\n",
    "import plotly.graph_objects as go\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "from scipy.optimize import minimize"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0119474c",
   "metadata": {},
   "source": [
    "## Import data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 397,
   "id": "f0e66963",
   "metadata": {},
   "outputs": [],
   "source": [
    "rho = np.linspace(0.1,6.1,31)\n",
    "Ls= np.arange(8,16,2).astype(float)\n",
    "\n",
    "Y = np.zeros((len(Ls),len(rho)))\n",
    "dY = np.zeros((len(Ls),len(rho)))\n",
    "for index, L in enumerate(Ls):\n",
    "    data = np.load('data/cluster/'+'2nn-{}-combined.npy'.format(int(L)))\n",
    "    Y[index] = np.mean(data,axis=0)\n",
    "    dY[index] = np.std(data,axis=0)\n",
    "    \n",
    "\n",
    "rho = rho[5:]\n",
    "Y = Y[:,5:]\n",
    "dY = dY[:,5:]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2e149d6",
   "metadata": {},
   "source": [
    "## scale data func."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 398,
   "id": "48745409",
   "metadata": {},
   "outputs": [],
   "source": [
    "def x_dat_2_new(rho, rho_c, l, nu):\n",
    "    return np.power(l, 1./nu) * (rho-rho_c)\n",
    "\n",
    "def x_BKT(rho, rho_c, l, nu):\n",
    "    return np.power(l* np.exp(abs((rho-rho_c))), 1./nu)\n",
    "\n",
    "def x_BKT(rho, rho_c, l, nu):\n",
    "    return np.power(l* np.exp(np.sqrt(abs(rho-rho_c))), 1./nu)\n",
    "\n",
    "def x_BKT(rho, rho_c, l, nu):\n",
    "    a = nu\n",
    "    return l*np.exp(-a/np.sqrt(abs(rho-rho_c)))\n",
    "\n",
    "\n",
    "def y_dat_new(y, l, zeta,nu):\n",
    "    \n",
    "    \n",
    "    return np.power(l, -zeta/nu) * y\n",
    "\n",
    "\n",
    "l_mesh_local, rho_mesh_local = np.meshgrid(Ls, rho, indexing='ij')\n",
    "\n",
    "def goal_function(rho_c = 3., nu=2.2, zeta=1.2):\n",
    "\n",
    "    #XX = x_dat_2_new(rho=rho_mesh_local, rho_c=rho_c, l=l_mesh_local ,nu=nu) \n",
    "    XX = x_BKT(rho=rho_mesh_local, rho_c=rho_c, l=l_mesh_local ,nu=nu) \n",
    "    \n",
    "    \n",
    "    YY = y_dat_new(y=Y, l=l_mesh_local, zeta=zeta, nu=nu)\n",
    "    #print('y',YY)\n",
    "    dYY = y_dat_new(y=dY, l=l_mesh_local, zeta=zeta, nu=nu)\n",
    "    #print('dy',dYY)\n",
    "    \n",
    "    # sorting for BKT\n",
    "    sorting_index = np.argsort(XX)[0]\n",
    "    #print(sorting_index)\n",
    "    \n",
    "    XX = XX[:,sorting_index]\n",
    "    YY = YY[:,sorting_index]\n",
    "    dYY = dYY[:, sorting_index]\n",
    "    \n",
    "    \n",
    "    return quality(XX,YY,dYY)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 399,
   "id": "4e932cca",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.2955174955843305"
      ]
     },
     "execution_count": 399,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "goal_function()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 400,
   "id": "00af36cc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "1.5\n",
      "2\n",
      "nu = [1, 1.5, 2] (vertical)\n",
      "zeta = [1.5, 2, 2.5] (horizontal)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAEYCAYAAABBWFftAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtj0lEQVR4nO3de5BkZ3nf8d9vdQEBzmpgJ1il1bC4iiIl8HDJcAsOkXHhCJCF03ZSkgOJ0OI1sWQuRrERKcPacdkEb2QuIuC1tBayZQSGwRZEBHAZLLC5zQoxIAkcQdAw8hqNpNFIWGstYp/88XZremfn0mfm9Hsu/f1UTfXtTPczXT1Pn+e9OiIEAAAAANi6bVUHAAAAAABtQYEFAAAAACWhwAIAAACAklBgAQAAAEBJKLAAAAAAoCQnVvXCO3bsiF27dlX18gCKOnxYWlyUjhyRTj5ZGhuTTjllw187ePDgXRExniHCdZFzgIbZRM4h3wDYlJLPcSorsHbt2qWZmZmqXh5AEbOz0r59KeFs3y4tLaVEdMkl0uTkur9q+/ZMUa6LnAM0yCZzDvkGQGFDOMdhiCCAjU1Pp8QzNiZt27Z8fXq6knBsP9L2l2x/1fbNtn+zkkAADEfNcg6AFhtCvimlwLJ9hu1P276le7Lz2jKeF0BNzM2lVp1+27en+6vxoKQXRsTTJD1d0tm2n1tVMABKVr+cI9sn2P6K7Y9VFgSA8g0h35TVg/WQpDdExJmSnivpIttnlvTcAKo2MZG6zPstLaX7KxDJ97s3T+r+RCXBAChfzXJO12sl3VplAACGYAj5ppQCKyIORcSN3ev3KyWg08t4bgA10Omk8ciLi9LRo8vXO53KQuq2Jt8k6U5Jn4qIL65yzB7bM7ZnFhYWsscIYJNqlnNs75T0UklXVBIAgOEZQr4pfQ6W7V2SniHpuJMdAA01OZkme46NSfPz6XKAyZ/DFBE/jIinS9op6dm2n7rKMfsjYioipsbHK19YDMCg6pdz3i7p1yQdXesAGnSAhhpCvil1FUHbj5H0YUmvi4j7Vnl8j6Q9kjRRbTc/gKImJystqNYSEffa/rSksyV9vep4AJSkJjnH9jmS7oyIg7bPWuu4iNgvab8kTU1NMWQZaJKS801pPVi2T1Iqrq6JiFWX3aA1GUAZbI/bPrV7/RRJL5L0jUqDAtBWz5d0ru3vSLpW0gtt/0m1IQGos7JWEbSkKyXdGhGXlfGcALCO0yR92vaspC8rzcFiZS8ApYuISyNiZ0TsknSepL+KiJdXHBaAGitriODzJb1C0te6k84l6U0RcX1Jzw8AD4uIWaW5ngAAoM1mZ9OeVHNzaWW/TqcWw4fXU0qBFRGfk+QyngsAAKCOIuIzkj5TcRjA6JidlfbtSwtP7NyZVvfbt6/qRW82VOoiFwAAoOUa2JoMoKGmp1NxNTaWbvcup6drnXdKX6YdAAC0VK81eXHx2Nbk2dmqIwPQRnNz0vbtx963fXu6v8YosAAAwGD6W5O3bVu+Pr3q4sEAsDUTE9LS0rH3LS2l+2uMIYJAGzBkB0AOc3Op56pfA1qTgVrhO3twnU7qJZdSrllaSj3nu3dXG9cGKLCApmvoBFCgNjjZGdzERMoxvXkQUiNak4HayPWdnSuvDft1JifTe9P/Grt31z5HM0QQaDqG7ACbx5yiYjqd9B4tLkpHjy5f73Sqjgxohhzf2bnyWq7XmZyU9u6VDhxIlzUvriQKLNTV7Gz6J7rwwnTJyc7aGjoBFKgFGiiK6bUmj41J8/Ppkt5yYHA5vrNz5TXy55oYIoj6YchbMQzZATaPOUXFTU6Si4HNyvGdnSuvkT/XRA8W6ocWkWIYsgNsXs4VqnL0zNP7D9Rbju/sXHmtoSv85UCBhfphyFsxDNkBNi9XA0WOuQrMJwPqL8d3dq68RgPvmhgiiPphyFtxDNkBNifXClX9PfPS8uX0dHmvleM1AGzdsL+zc+W1hq7wlwMFFuqnoXseIB/bZ0i6WtLjJYWk/RHxjmqjQmPlaKDIMVeB+RAAenI1vNLAuyqGCKJ+GPKGjT0k6Q0Rcaak50q6yPaZFccErC3HXAXmQwBALdCDhXqiRQTriIhDkg51r99v+1ZJp0u6pdLAgLXk6Jmn9x8AaoEeLGDYWNVrqGzvkvQMSV9c5bE9tmdszywsLGSPDXhYjp55ev8BoBbowQKGiT29hsr2YyR9WNLrIuK+lY9HxH5J+yVpamoqMocHHCtHzzy9/wBQOQosjK7Z2WNXvul0mrly2IiyfZJScXVNRLBJGgBg63KcG6D1GCKI0ZRrvxj29BoK25Z0paRbI+KyquPBEDHEFkAu7CWHklBgYTT19yxt27Z8fbrkjhBW9RqW50t6haQX2r6p+/OSqoNCyTjZAZBTrnMDtB5DBFFcG7rPc+0Xw6peQxERn5PkquPAkDHEFkBO7CWHktCDhWLa0qKcq2eJVb2AzWOILYCcGHWCktCDhWLa0qKcs2eJVb2AzZmYSP+XvTwjcbIDYHgYdYKS0IPVJjkmg7elRZmeJaD+Op10crO4KB09uny906k6MgBtxLkBSlJaD5btA5LOkXRnRDy1rOfFgHLtt9SmFmV6loB6653s9M/53L2b/1tkZfsMSVdLerykkLQ/It5RbVQYGs4NUIIyhwheJelypSSE3HIN3aP7HEBOnOygeg9JekNE3Gj7RyQdtP2piLil6sBGShsW2MLIKG2IYETcIOmesp4PBeUaukf3OYAe9qjCCIiIQxFxY/f6/ZJulXR6tVGNmLYssIWRkXWRC9t7JO2RpIkmDimrs5xD92hRBpBrWDJQI7Z3SXqGpC+u8hjnOMPSlgW2MDKyLnIREfsjYioipsbHx3O+dPsxGRxATmzIiRFj+zGSPizpdRFx38rHOccZorYssIWRwSqCbcHQPQA5ccKDEWL7JKXi6pqIoBUhN/anQsOwD1abMHQPQC5tWlEUWIdtS7pS0q0RcVnV8YwkFthCw5TWg2X7/ZI+L+nJtudt86kHgLZiWDJGx/MlvULSC23f1P15SdVBjRRG6aBhSuvBiojzy3ouAEDNsUcVRkREfE6Sq45j5DFKBw3CEEEAjcPG5jXBCQ8AAMdhkQsATXSVpLOrDgIAAGAlCiwAjcPG5gAAoK4osAC0lu09tmdszywsLFQdDgAAGAEUWABai40/AQBAbhRYAAAAAFASVhEEAADA5s3OHrtlQ6fDCqMYafRgAWgcNjYHgJqYnZX27Usbje/cmS737Uv3AyOKHiwAjcPG5gOgRRlADtPT0thY+pGWL6enyTkYWRRYuXCyAyCXXovy2NixLcqXXELeAVCuubmUZ/pt357uB0YUQwRzoPscQE79Lcrbti1fn56uOjIAbTMxIS0tHXvf0lK6HxhRFFg5cLIDIKe5udSC3I8WZQDD0OmkhuPFReno0eXrnU7VkQGVocDKgZMdADnRogwgl8nJNPx4bEyan0+XDEfGiGMOVg4TE6k1pzfxU+JkB8DwdDppGLKUGnOWllIO2s1iiwCGYHKSggroU+8Cqy0LQ3CyAyCnXotyf/7cvbuZ+RPA1rTlXApokPoWWG1aBYuTHQC50aIMoE3nUkCD1LfAatu+CpzsAJBoTQaQT9vOpYCGqO8iFywMAaBt2LIBQE6cSwGVqG8PVs6FIWhRBiANPxfQmgygJ8e5B4tsAZWobw9Wrn0VaFEGIOXJBbQmA5DynXuwRxVQifoWWLn2VWATYABSnlzA/lQApHznHuxRBVSivkMEpTwLQ8zNpdajfrQoA6MnRy5gywYAUt5zDxbZArKrbw9WLrQoA5Dy5AJakwFInHsALVfvHqwcaFEGIOXLBbQmA+DcA2g1R0Q1L2wvSLp9wMN3SLprWLE8SjrlsdKpJ0uPOCI9eI907wPS4arj2gLiKoa4iika1xMiYnxYwQxqkJzTywX/JG1/pLRUIBfk0pbPRC7EVUwb4iLflKcNn4eciKuYtsS1as6prMAqwvZMRExVHcdKxFUMcRVDXNWp699IXMUQVzHEVY26/n3EVQxxFdP2uJiDBQAAAAAlocACAAAAgJI0pcDaX3UAayCuYoirGOKqTl3/RuIqhriKIa5q1PXvI65iiKuYVsfViDlYAAAAANAETenBAgAAAIDao8ACAAAAgJJQYAEAAABASSiwAAAAAKAkFFgAAAAAUBIKLAAAAAAoCQUWAAAAAJSEAgsAAAAASnJiVS+8Y8eO2LVrV1UvD6Cow4elxUXpyBHp5JOlsTHplFM2/LWDBw/eFRHjGSJcFzkHaJhN5BzyDYBNKfkcp7ICa9euXZqZmanq5QEUMTsr7duXEs727dLSUkpEl1wiTU6u+6u2b88U5brIOUCDbDLnkG8AFDaEcxyGCALY2PR0SjxjY9K2bcvXp6erjgxAG5FzAOQyhHxDgQVgY3NzqVWn3/bt6X4AKBs5B0AuQ8g3FFgANjYxkbrM+y0tpfsBoGzkHAC5DCHfUGAB2Fink8YjLy5KR48uX+90qo4MQBuRcwDkMoR8Q4EFYGOTk2my59iYND+fLgeY/AkAm0LOAZDLEPJNZasIAmiYyUlObgDkQ84BkEvJ+YYeLAAAAAAoCQUWAAAAAJSEAgsAAAAASkKBBQAAsAbbZ9j+tO1bbN9s+7VVxwSg3iiwANSK7QO277T99TUeP8v2ku2buj9vzh0jgJHykKQ3RMSZkp4r6SLbZ1YcE4AaYxVBAHVzlaTLJV29zjGfjYhz8oQDYJRFxCFJh7rX77d9q6TTJd1SaWDAqJidlaanpbm5tPlvp1P7FUbpwQJQKxFxg6R7qo4DAFayvUvSMyR9cZXH9tiesT2zsLCQPTaglWZnpX370sa/O3emy3370v01RoEFoImeZ/urtj9u+ylrHcQJD4Cy2H6MpA9Lel1E3Lfy8YjYHxFTETE1Pj6eP0Cgjaan08a/Y2PStm3L16enq45sXRRYAJrmRklPiIinSXqXpD9f60BOeACUwfZJSsXVNRFR7zM7oE3m5qTt24+9b/v2dH+NMQcLQKP0txxHxPW2/5ftHRFxV5VxASOjgfMhtsK2JV0p6daIuKzqeNACI/Y/tCUTE2lY4NjY8n1LS+n+GqMHq01mZ6W9e6ULL0yXNR+fCmyG7R/tnvDI9rOV8tjd1UYFjIiGzofYoudLeoWkF/atXvqSqoNCQ43m/9DmdTrpPVpclI4eXb7e6VQd2bo27MGyfUDSOZLujIinrvK4Jb1D0kskPSDpgoi4sexAsYHeP+zY2LH/sJdcQqsIGsX2+yWdJWmH7XlJb5F0kiRFxHsl/byk/2L7IUmHJZ0XEVFRuGgDWpMH1z8fQlq+nJ5u7XsWEZ+T5KrjQEuM4P/QlkxOpnPZ/hy9e3ft36tBhghepfWXTH6xpCd1f54j6T3dS+TEPyxaIiLO3+Dxy5VyErB1NE4VMzeX3qd+DZgPAdQG/0PFTU42Lh9vOERwgCWTXybp6ki+IOlU26eVFSAG1NBJgABQqYauUFWZiYk0/6FfA+ZDALXB/9BIKGMO1umSvtt3e75733FYMnmI+IcFgOJonCqmofMhgNrgf2gkZF3kgiWTh6ht/7As2AEgBxqniunNhxgbk+bn0yXDKYHBte1/iPO1VZWxTPsdks7ou72zex9yaugkwFUxJwJALp1Oyi9S6rlaWko5Z/fuauOqswbOhwBqJcf/UI7FezhfW1MZBdZ1ki62fa3S4hZLEXGohOdFUW350mPBjuJYBQ3YnDY1TgGAlK/w4XxtTYMs077RksnXKy3RfpvSMu2vHFawGBGssFMMLUjA1rSlcQoApHyFD+dra9qwwBpgyeSQdFFpEQEN3bW7MrQgAQCAnlyFD+dra8q6yAUwkLYt2DFsrIIGAAB6ci3ew/namsqYg4VRM+z5PsyJKIYWJAA9zMcEkGvxHs7X1kSBhWJyzfdhTsTgWAUNbUbBMDjmYwJb14ack7Pw4XxtVRRYKIb5PvXTshYk2wcknSPpzoh46iqPW9I7lBbXeUDSBRFxY94oIWn4JyJtKxiG/X6Rn4GtaVPOofCpFHOwUAzzfeppcjJt8HfgQLpsdlK9StLZ6zz+YklP6v7skfSeDDFhpd6JyOLisSciZW4y2V8wbNu2fH16urzXyCXH+0V+BramTTkHlaLAQjG5Jk5iZEXEDZLuWeeQl0m6OpIvSDrV9ml5osPDcpyItKlgyPF+kZ+BrWlTzkGlKLBQDCvGoHqnS/pu3+357n3Hsb3H9oztmYWFhSzBjYwcJyJtKhhyvF/kZ2Br2pRzUCkKLBTTm+8zNibNz6fLJo5NxkiIiP0RMRURU+Pj41WH0y45TkTaVDDkeL/Iz8DWtCnnoFIscoHimDiJat0h6Yy+2zu79yGnHKtXtmkBl5zLJjfx/QE2kmN1vzblHFSKAgtA01wn6WLb10p6jqSliDhUcUyjJ9eJSFsKBk7cgM3LubpfW3IOKkWBBaBWbL9f0lmSdtiel/QWSSdJUkS8V9L1Sku036a0TPsrq4kUnIgUxPvVWBttH4EhYwsCNAwFFoBaiYjzN3g8JF2UKRwAkNL2EZdLurriOEbT3FzquerH6n6oMRa5AAAAWMcA20dgmFjdDw1DgQUAALBFbAsxRKzuh4ahwAIAANgitoUYIrYgQMMwBwsYthxLywIA0GYsEoMGqXeBxYkpmi7n0rIAAACoXH2HCPZOTBcXjz0xnZ2tOjJgcP1Ly27btnx9errqyAAAA+puH/F5SU+2PW+75B2iAbRJfXuw2PMAbcDSsgDQeBttHwEA/erbgzU3l05E+3FiiqZhaVkAAICRUt8erImJNCyw13MlcWKKcuWY49fppKGtUmogWFpKn+vdjC4BAABoo/r2YLHnAYYp1xw/lpYFAAAYKfXtweqdmPb3MOzezYkpypFzjh9LyxZi+2xJ75B0gqQrIuKtKx6/QNLvSbqje9flEXFF1iABAADWUN8CS2rXiSlLztcLi0/Uku0TJL1b0oskzUv6su3rIuKWFYd+ICIuzh4gAADABuo7RLBNWHK+flh8oq6eLem2iPh2RByRdK2kl1UcUzPNzkp790oXXpguyTcAAGQxUIFl+2zb37R9m+03rvL4BbYXbN/U/XlV+aE2GHsh1Q9z/OrqdEnf7bs9371vpZ+zPWv7Q7bPWOvJbO+xPWN7ZmFhoexY64tGHQA50aADHGPDAqtvyM6LJZ0p6XzbZ65y6Aci4undH+ZD9GPJ+fph8Ykm+6ikXRExKelTkt631oERsT8ipiJianx8PFuAlaNRB0AuNOgAxxlkDtbDQ3YkyXZvyM7KORFYC0vO11Ob5vi1xx2S+nukdmp5MQtJUkTc3XfzCklvK+3V2zJXkjmGAHLJuWgU0BCDDBEsbcjOyA7XYTgaMKgvS3qS7SfaPlnSeZKu6z/A9ml9N8+VdGspr9ymVljmGALIhVE6wHHKWuRioCE7Iztch+FowEAi4iFJF0v6hFLh9MGIuNn2b9k+t3vYa2zfbPurkl4j6YJSXrxNw+po1AGQCw06wHEGGSJY7ZCdtmA4GjCQiLhe0vUr7ntz3/VLJV1a+gu3aVgd+wgCyKXTSb39UsqZS0upQWf37mrjAio0SIH18JAdpcLqPEm/0H+A7dMi4lD3ZnlDdgAgl7bNlaRRB0AONOgAx9mwwIqIh2z3huycIOlAb8iOpJmIuE5pyM65kh6SdI/KGrIDALnQCgsgp7YsqiPRoAOsMEgPVnVDdnJpU5IDsDm5WmHJNwB6i+qMjR27qM4w5meTc4DsylrkornatHIYgK2ZnEybZB44kC6HcaJDvgGQa1Edcg5QCQqsNq0cBqDeyDcApHxLm5NzgEpQYLF/A4BcyDcApHxLm5NzgEpQYLF/A4BcyDcApHx71ZFzgEpQYLEhJ4BcyDcApOVFdcbGpPn5dDmMBS7IOUAlHBHVvLC9IOn2AQ/fIemuYcXyKOmUx0qnniw94oj04D3SvQ9Ih6uOawuIqxjiKqZoXE+IiPFhBTOouuScLeSboca1RcRVDHEVUyQu8s0KnONkQ1zFtCWuVXNOZQVWEbZnImKq6jhWIq5iiKsY4qpOXf9G4iqGuIohrmrU9e8jrmKIq5i2x8UQQQAAAAAoCQUWAAAAAJSkKQXW/qoDWANxFUNcxRBXder6NxJXMcRVDHFVo65/H3EVQ1zFtDquRszBAgAAAIAmaEoPFgAAAADUHgUWAAAAAJSEAgsAAAAASkKBBQAAAAAlocACAAAAgJJQYAEAAABASSiwAAAAAKAkFFgAAAAAUBIKLAAAAAAoyYlVvfCOHTti165dVb18Ox0+LC0uSkeOSCefLI2NSaecUnVUGHEHDx68KyLGq46DnAM0zCa+08g3ADZlk+fQa+WcygqsXbt2aWZmpqqXb5/ZWWnfvvSB2L5dWlpKH5RLLpEmJ6uODiPM9u1VxyCRc4BG2eR3GvkGQGFbOIdeK+cwRLAtpqfTB2NsTNq2bfn69HTVkQEAUAzfaQByGUK+ocBqi7m5VHX327493Q8AQJPwnQYglyHkGwqstpiYSF2a/ZaW0v0AADQJ32kAchlCvqHAaotOJ40XXVyUjh5dvt7pVB0ZAADF8J0GIJch5BsKrLaYnEyT8cbGpPn5dMkCFwCAJuI7DUAuQ8g3la0iiCGYnOTLBwDQDhm+02wfkHSOpDsj4qmrPH6WpL+Q9P+6d01HxG8NNSgA+ZWcbyiwAADAqLpK0uWSrl7nmM9GxDl5wgHQBgwRBAAAIykibpB0T9VxAGiXDQss2wds32n762s8btvvtH2b7Vnbzyw/TAAAgEo8z/ZXbX/c9lPWOsj2HtsztmcWFhZyxgegZgbpwbpK0tnrPP5iSU/q/uyR9J6thwUAAFC5GyU9ISKeJuldkv58rQMjYn9ETEXE1Pj4eK74ANTQhgXWAN3nL5N0dSRfkHSq7dPKChAAAKAKEXFfRHy/e/16SSfZ3lFxWABqroxFLk6X9N2+2/Pd+w6tPND2HqVeLk2wWSDQLLOz0vR02tl8YiLtD8GqlQBazPaPSvpeRITtZys1TN9dcVgAai7rIhd0nwMNNTsr7duXNt7buTNd7tuX7geAhrL9fkmfl/Rk2/O2d9t+te1Xdw/5eUlft/1VSe+UdF5ERFXxAmiGMnqw7pB0Rt/tnd37ALTF9HTaeG9sLN3uXU5P04sFoLEi4vwNHr9caRl3ABhYGQXWdZIutn2tpOdIWoqI44YHbgpDkoB6mJtLPVf9tm9P9wMAAOBhGxZY3e7zsyTtsD0v6S2STpKkiHivpOslvUTSbZIekPTKUiLrDUkaGzt2SNIll1BkAblNTKT/wV7PlSQtLaX7AQAAhqWBHS4bFlgDdJ+HpItKi6iHIUlAfXQ6qYFDSj1XS0up4Nq9u9q4UK0GfukBABqkoR0uWRe5KGRuLp3I9WNIElCNycmUzMbGpPn5dDmk5Mbm5g3BwicAsDmzs9LevdKFF6ZL8uba+jtctm1bvj49XXVk6ypjDtZwMCQJqJfJyVytRVcpTSq/eo3H+zc3f47S5ubPyREY+jDKAACKa2iPTGUaOge8vj1YnU760C0uSkePLl/vdKqODMAQsbl5QzDKAACKa2iPTGUmJlIHS78GdLjUt8DKOCQJQKOstbn5cWzvsT1je2ZhYSFLcCOjoV96q2K4DoBcaJwqpqEdLvUdIijlHJIEoIUiYr+k/ZI0NTXF5qBlasvCJwzXAZATU2CK6XW49C+otHt37fNzvQssADgem5vXQUO/9I7DXLLiWD0S2Ly2NE7l1MAOFwosFMeXK6o1vM3NUUwDv/SO09AJ1JWhxw/YmrY0TmFdFFgohi9XDFllm5tjNDFcpxh6/ICta0PjFNZFgYVi+HLFkFW2uTlGE8N1iqHHDwA2VN9VBFFPrH4DoE1YsbaYNq0eKTY2BzAc9GChmFzDaZjnBUDKkwsYrjO49vX4XSU2Nm8GzgvQIPRgoZgc+xH05nktLh47z4u9aYDRQi6on5b1+LGxeUPkzAXsi4cS0IOFYnKsfsM8LwASuaCuRqvHb62NzY9budT2Hkl7JGmioUMmaytXLmAhr+LoWVwVBRaKG/aXK5OoAUjkAjQKG5sPUa5cQKNOMRSka2KIIOqnZZOoAWwSuQDVY2PzOsiVC1jIq5j+gnTbtuXr09NVR1Y5CizUT455XgDqj1yA6l0n6T91VxN8rtjYvBq5ckGbGnVyzCWjIF0TBVYuTJocXMsmUQPYJHIBhqy7sfnnJT3Z9rzt3bZfbfvV3UOul/RtpY3N/1DSL1cU6mjLlQva0qiTa1GQNhWkJWMOVg6MUS1utCZRA1gLuQBDxMbmDZIjF+RYyCuHXHPJ2rdtQ2kosHJg0iSGjVV8AADYujY06uRaFKQtBekQUGDlwEpYGCZ6SAEAQM/ERDoX6DXoS8MbuteGgnQImIOVA2NUMUys4gMAAHraMpeswSiwcuCDjmFiFR8AANDDAkGVY4hgDoxRxTDlHAoAAADqj6F7laLAkvIsEMAHHcPCKj4AAAC1QYHFAgFoOnpIy8NqjGg6PsMAUDkKLJZQRxvQQ7p1NLag6fgMA0AtsMgFCwQAkFiNEc3HZxgAaoECiyXUAUg0tqD5+AwDQC1QYLGEOgCJxhY0H59hAKgFCiz2CgAg0diC4ZudlfbulS68MF3Ozpb7/HyGAaAWBiqwbJ9t+5u2b7P9xlUev8D2gu2buj+vKj/UIZqcTF92Bw6kS4orYPTQ2IJh6i1Asbh47AIUZRZZfIY3pfXnOACy23AVQdsnSHq3pBdJmpf0ZdvXRcQtKw79QERcPIQYASAPVmPEsORasZbPcCGc4wAYhkGWaX+2pNsi4tuSZPtaSS+TtDL5AADqgv2Q6mVuLvVc9WMBijrgHAdA6QYZIni6pO/23Z7v3rfSz9metf0h22es9kS299iesT2zsLCwiXABtB3DdUqQYzgaimEBiroq7RxnpA17fiHQMGUtcvFRSbsiYlLSpyS9b7WDImJ/RExFxNT4+HhJLw2gLfqG67xY0pmSzrd95iqHfiAint79uSJrkE3Afkj1wwIUTTbQOU7hRuS2FCU06ADHGaTAukNSf2vNzu59D4uIuyPiwe7NKyT9y3LCAzBiHh6uExFHJPWG66AI9kOqHxagqKvSznEKNSK3qSihQQc4ziBzsL4s6Um2n6iUdM6T9Av9B9g+LSIOdW+eK+nWUqMEmoy5MEWsNlznOasc93O2XyDp7yS9PiK+u8oxo2tiIp2w9RZSkBiOVgcsQFFH1Zzj5Fr0JAfmFwLH2bAHKyIeknSxpE8oJZUPRsTNtn/L9rndw15j+2bbX5X0GkkXDCtgoFHa1EpZHwMN15FGeN4nw9GAgVR2jtOmXmbmFwLHGaQHSxFxvaTrV9z35r7rl0q6tNzQgBZoUytlHgMN1+m7eYWkt631ZBGxX9J+SZqamorywqy53nC0/p7T3bv5zAGrqOQcp029zJ1OajiUUpG4tJT+tt27q40LqNBABRaATWLoRFEMSS4Lw9GA+mpTUUKDDnAcCixgmNrUSplBRDxkuzdc5wRJB3rDdSTNRMR1SsN1zpX0kKR7xJBkAE3TtqKEBh3gGBRYwDC1qZUyE4YkAxgJFCVAa1FgYXTlWN2vba2U2DpWlQSQEzkHyI4CC6Opt7rf2Nixq/sNY18aWinRk/NzBwDkHKASg2w0DLQPGyOiCnzuAOREzgEqQYGF0dSmPUjQHHzuAOREzgEqQYGF0cTGiKgCnzsAOZFzgEpQYGE0dTppLPrionT06PL1TqfqyNBmfO4A5ETOASrhiKjmhe0FSbcPePgOSXcNMZzNIq5iahXXo6RTHiud+k/S9kdKS/dI9z4gHa46rj61er/6FI3rCRExPqxgBlWXnNP73J0sPeKI9GDBz11bPhO5EFcxbYiLfLPCFnJOGz4PORFXMW2Ja9WcU1mBVYTtmYiYqjqOlYirGOIqhriqU9e/kbiKIa5iiKsadf37iKsY4iqm7XExRBAAAAAASkKBBQAAAAAlaUqBtb/qANZAXMUQVzHEVZ26/o3EVQxxFUNc1ajr30dcxRBXMa2OqxFzsAAAAACgCZrSgwUAAAAAtUeBBQAAAAAlqXWBZfuA7Tttf73qWPrZPsP2p23fYvtm26+tOiZJsv1I21+y/dVuXL9ZdUz9bJ9g+yu2P1Z1LD22v2P7a7Zvsj1TdTw9tk+1/SHb37B9q+3n1SCmJ3ffp97PfbZfV3VcZapjziHfbA75ZnDkm2rUMd9I5JzNIucMbhRyTq3nYNl+gaTvS7o6Ip5adTw9tk+TdFpE3Gj7RyQdlPSzEXFLxXFZ0qMj4vu2T5L0OUmvjYgvVBlXj+1flTQl6Z9FxDlVxyOl5CNpKiJqtdmd7fdJ+mxEXGH7ZEmPioh7Kw7rYbZPkHSHpOdExKCbadZeHXMO+WZzyDeDI99Uo475RiLnbBY5Z3CjkHNq3YMVETdIuqfqOFaKiEMRcWP3+v2SbpV0erVRSZF8v3vzpO5PLSpo2zslvVTSFVXHUne2t0t6gaQrJSkijtQp8XT9lKRvtelkR6pnziHfFEe+GRz5pjp1zDcSOWczyDmDG5WcU+sCqwls75L0DElfrDgUSQ93Ud8k6U5Jn4qIWsQl6e2Sfk3S0YrjWCkkfdL2Qdt7qg6m64mSFiT9UXe4wRW2H111UCucJ+n9VQcxasg3A3u7yDeDIt9gTeScgb1d5JxBjUTOocDaAtuPkfRhSa+LiPuqjkeSIuKHEfF0STslPdt25cMObJ8j6c6IOFh1LKv4iYh4pqQXS7qoO2SjaidKeqak90TEMyT9o6Q3VhvSsm53/rmS/qzqWEYJ+WYw5JvCyDdYFTlnMOScwkYi51BgbVJ3/O+HJV0TEdNVx7NSt7v105LOrjgUSXq+pHO7Y4GvlfRC239SbUhJRNzRvbxT0kckPbvaiCRJ85Lm+1rmPqSUjOrixZJujIjvVR3IqCDfFEK+KYZ8g+OQcwoh5xQzEjmHAmsTuhMtr5R0a0RcVnU8PbbHbZ/avX6KpBdJ+kalQUmKiEsjYmdE7FLqdv2riHh5xWHJ9qO7E3jV7Z7+aUmVr+YUEf8g6bu2n9y966ckVTq5eIXzxXCdbMg3xZBviiHfYCVyTjHknGJGJeecWEIgQ2P7/ZLOkrTD9rykt0TEldVGJSm1VrxC0te6Y4El6U0RcX11IUmSTpP0vu7qJ9skfTAiarNcaA09XtJH0neJTpT0pxHxf6oN6WG/Iumablf1tyW9suJ4JD2cpF8k6ZeqjmUYappzyDftQL4piHxTGXJOO5BzCioz59R6mXYAAAAAaBKGCAIAAABASSiwAAAAAKAkFFgAAAAAUBIKLAAAAAAoCQUWAAAAAJSEAgsAWs7282x/0Pbf2z5i+27bn7L9n7tLHm/0+1d1N9Hs3d5le6/tH1vl2O/Yvqrcv+Dh5/6M7c8McFz0/fzA9oLtz9r+Ddv/fLPPm0M35r1VxwEA2Lxa74MFANga26+TdJmkv5L065JulzSmtOnkeyTdK+kvNnia/y7pHX23d0l6i6TPKe1h0u/fSbpva1GX4ipJf6DUkPg4Sc9V2nvlNbZfFhF/23fsL+cPDwDQVhRYANBStl+gVFxdHhGvWfHwX9i+TNKj1/n9R0TEgxHxrUFfMyK+srloS3dHRHyh7/ZHbb9T0mclTdv+sYh4QJIi4pZcQfXe07a9FgBgGUMEAaC9fl3SPZJ+bbUHI+JbETErSbYv6A5Pe4HtP7N9r6Qvdh97eIig7bMkfbr7FJ/qG4p3Vvfx44YI2n6i7T+2/Q+2H7T9bdvv6Hv8WbY/ZHve9mHb37T9O7ZPKeuN6P6935P0XyU9XtL5fa9/zBBB24+x/S7bc91477T9l7b/Rd8x/8z25d1hlw92Y369bfcdc1b3venY/kPbC5K+133sBNu/bfuQ7Qe6MTxltbhtP832dbYXu+/P39j+1yuOuar7/j3P9t/aPizpbeW8cwCAIujBAoAW6s6t+klJfx4R/1TgV6+R9H5JP6/VvyNulHSRpHdLeo2kL3fvX7UXyPYTJX1J0gOS3izp/0qaUBqi2DMh6SalYX33S3pK99gfk3RegdgH8UlJD0l6vqQr1zjm9yWdK+lN3Xgf1z3+VEmyvU3S/5b0zG6cX5P0UqXewvHu7/V7l6SPS3qFpEd279vbPe6ybkxTkq5bGYjtZyr1un1F0i8qvY+vlvSXtv9VRBzsO3y7pGsl7es+9+F13wkAwFBQYAFAO+2QdIrSnKsiPhQRq/Z4SVJE3Ge7V0zdumIY3mp+sxvH0yLi7/vuf1/fc364d73bA/Q3SvO4rrZ9UUTcXfBvWFNEHLZ9l6TT1jnseZKuiYj+AuwjfddfIuknJL0yIq7q3vdJ24+W9Abbl0XEXX3HfykiXtW7YXtM0usl7Y+IS/p+/4eS3roilt+TNCfphRFxpPv7n5D0dUm/Ieln+459jKSXR8RGc+oAAEPEEEEAQL+PbHxIIT8t6WMriqtjdIfb/Q/b35L0oKQfSPpjSZb0pJLjUfd5Y53HvyzpAttvsj21ykqLL5B0VNKfrrj/TySdrFSg9Vv5nv640ty3D664/9pjgkxDJP+NpD+TdNT2ibZP7Mb/l904+v1A0sfW+bsAABlQYAFAO92tNETsCQV/71DJcTxO0vwGx/yR0rC3d0p6kaRnKQ1DlJaH1JWiW7Ts0Pp/568orUB4oVKxdaft37f9qO7jj5V0T69Hqc8/9D3eb+Vr9XrPvrfi/pW3HyvpBKWeqh+s+LlY0lh3uGLPQkT8cJ2/CwCQAUMEAaCFIuKh7sINLyq4mtx6PTubcZek09d60PYjJb1M0t6I6F/44sdLjqPn3yoVLZ9b64CI+L6kSyVdavsJSvPR3irpiJYXDnms7ZNXFFk/2r28Z+VTrrjdK7geL+nmvvsfv+K4e5V6yt4t6eo1Yj26zusAACpADxYAtNdblXqQVl1Nrru63+QmnrdXrA2yyt8nJZ1je605T49QKnh+sOL+CzYR17q6mwy/TanAuXaDwyVJEXF7RPxPpYUsntq9+6+Vvj///YrD/6NSEfb5DZ52VtI/SvoPK+4/ZkGPiPhHpQUunibpxoiYWfkzyN8AAMiLHiwAaKmIuMH2r0q6zPaZSqv0zSltNPxTkl4l6ReUTviL+DullfgutH2PUsH1zYi4f5Vj36K0KMTf2v4dSbcp9WidHREvj4gl219QWhzikFKP14Vap9drQKfbfq5SIfRYpY2Gf1Fp/tLPRMSaK+zZ/rzSin5fk/R9pXlQT9PywhwfV+oBe6/tcaVeqJcovZ+/u2KBi+NExL22f1/Sf7N9v1IR+ixJu1c5/Fcl3SDpE7avVCoOdyitYHhCRLxxozcCAJAXBRYAtFhEvN32l5RWrdundHJ+v6QZSb8k6aObeM67bV+sNFzur5V6oH5S0mdWOfY73ULntyX9rtJKd3dI6l/p7nxJ71EaCndYafGH12prCzZc0P15SNKSpG8oLZf+BxGxsMHv3qDUu/RGpe/Jb0t6fUS8s/s3HbX9Ukm/o/QePE7Sd5SKobcPGN9epWLvVUrzqb4o6Wd07JBBRcSNtp+lVKi+U2kp9gWl5fLfO+BrAQAycgRDtgEAAACgDMzBAgAAAICSUGABAAAAQEkosAAAAACgJBRYAAAAAFASCiwAAAAAKAkFFgAAAACUhAILAAAAAEpCgQUAAAAAJfn/48OuInYbB4QAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x288 with 9 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(3,3, figsize=(12,4), sharex=True)\n",
    "\n",
    "for index0, nu in enumerate([1,1.5,2]):\n",
    "    print(nu)\n",
    "    for index1, zeta in enumerate([1.5,2,2.5]):\n",
    "        for r in np.linspace(1,7,10):\n",
    "            ax[index0,index1].scatter(r,goal_function(r,nu,zeta), c='r', alpha=.5)\n",
    "            \n",
    "print('nu = [1, 1.5, 2] (vertical)')\n",
    "print('zeta = [1.5, 2, 2.5] (horizontal)')\n",
    "\n",
    "ax[2,1].set_xlabel('Critical Disorder', fontsize=16)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "            \n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 487,
   "id": "5660eab2",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 10000/10000 [01:21<00:00, 122.30it/s]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fab2cdfc490>]"
      ]
     },
     "execution_count": 487,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAq2ElEQVR4nO3deZhV1ZU28HdVMckggwVRJgFFBURBSsB2QlFB1BinqEmMMXZIOklHzZektdtOUGN3tKNGOxokapxR40gcYjAOGBBlUEAElJkCZJ6Kmar9/bFq9zn33Hk49wz3/T1PPWe4t+5dBVXr7rPP3muLMQZERBR9VUEHQEREpcGETkQUE0zoREQxwYRORBQTTOhERDHRLKg3rqmpMb169Qrq7YmIImnWrFkbjTGdUz0WWELv1asXZs6cGdTbExFFkoisSPcYu1yIiGKCCZ2IKCaY0ImIYoIJnYgoJpjQiYhiggmdiCgmmNCJiGKCCZ2IqJw++wzYvt2Xl2ZCJyIqlzVrgOeeA154wZeXZ0InIiqXZct027+/Ly/PhE5EVC4bNgAtWwJDhvjy8kzoRETlsnQpcOSRQDN/ymgxoRMRlcP27frVs6dvb8GETkRUDmvW6LZbN9/eggmdiKgcVq8GqqqAQw/17S2Y0ImIymH1auArX/Gt/xxgQiciKo8NGzSh+4gJnYjIb42NQH09cPDBvr4NEzoRkd82bQKMATp18vVtmNCJiPy2ZIlufRzhAjChExH5769/1e0hh/j6NkzoRETlUuVvyvVv/AwRUaXbtQuYPVv3R43y/e2Y0ImI/PLYY8C6dbrv45R/K2v7X0QeEZH1IvJpmse/KSJzRWSeiEwTkeNLHyYRUcQsW+YkcwDo3Nn3t8ylQ+dRAKMzPL4MwOnGmIEAbgMwoQRxERFFx6ZNwMKFzvGBA8CzzyY+p0UL38PI2uVijJkiIr0yPD7NdTgdQPcSxEVEFB3/+7+6HTdOt7NnA3v2ABdcAKxfDxxzTFnCKHUf+rUA3ijxaxIRhZcxyedWrAA6dPBtIYt0SpbQReQMaEI/JcNzxgIYCwA9y3CDgIjId1u3Ovv79ul2/nzfJxGlUpJBkSJyHICHAFxojNmU7nnGmAnGmFpjTG3nMtwgICLyXV2ds//ss8CCBbrft2/ZQyk6oYtITwAvArjKGPN58SERhcDOncB776W+nCZyW7XK2V+yBHjpJd0fOrTsoWTtchGRiQBGAKgRkToAvwLQHACMMeMB/BLAIQAeEBEAOGCMqfUrYKKyePll4IsvgD59gB49go6GwqyuDujSRW9+WlddBbRuXfZQchnlcmWWx/8ZwD+XLCKiMPjiC90+/LAzcoHIq7FRx5rX1joJvUUL4IgjAgmHtVyIiApVXw80NAAdOzrnzjgjsHCY0Ily1dAQdAQUNl9+qduaGuece7/MmNCJUrEtroMO0u38+cBttwEbNwYXE4XP8uVAdXXifZY+fQILhwmdKJUjj9Tt8OG6fest3U6dGkw8FE7Ll2syb95cj9u31wQfECZ0okwaG3W7ZYtuP/44uFgoXPbuBdauBXr10jrn3/0u8MMfBhoSy+cSpWL7yxsbgR07nPM6NJdI5yoY43TPhWD2O1voRKnYhN7QANx1l3M+wMtpChn7OxKi3wkmdKJUbFeL3Vrt25c/FgonJnSiiHC30N1STedetAjYts3/mChcmNCJIiJdQvdavRqYOBG45x7/Y6JwYUInigj3TVE3b7Gu5cvLEg6F0M6dum3VKtg4XJjQqbIdOKDLh6U6D2RvoW/fXvqYKBps7ZYuXYKNw4UJnSrbK6/o8mF2YQIrXZeLMfrc/fv1uAzrRFJIrV2rqxLZSUUhwIROlc1WVfQmbntsW+qWMcB//Rdw++163LKlbo87zr8YKXx27tTfnUMPDTqSBEzoVNm8feRWuoRuW+be51XxT6mivPqqfrifeGLQkSTgbyFVtnQ3P+35tWsTz7/zTuJxug8Eiq/GRl1mbsCAwOqep8OETpXNJmTv6BXbMt+9u7Dvp/iaO1e3HToEGkYqrOVClc0m4jlzgJNPds7nWvucCb1y7Nun90+sABeySIctdCIgeaZnLgl92TIm9EqyerWzP2RIKO+bsIVOBACdOyce55LQH3sM6NpV95nQ4+3LL/X/2zr//OBiySB8HzFE5WRn+XXqlHg+1y6XNWt0y4Qeb+PHO/sDB4a2jDJb6FTZUiVuY5KHK2bDhB5fdtLZkUcCl18ONAtv2mQLnfxz1106XjfMUtVqmTw5/9dhQo8vO9Kpf3+dFRrS1jnAhE5+2rEDmDlT9+vrgcWLg40nFe9NzYYGYNo03ff2q2fChB5fe/boNkRFuNJhQqfy+O1vgSefDFfiMyY5Hvexe+m5XF6L4sneJ2nTJtg4csCETv5Il+BSnQ8qGW7cmBjDnj3Ahg3OucGDddu9e/bXYkKPp717tYBbTU0o1gzNhgmd/JFuBR/vzcYvvgBuuUWnUpeb94bon/4ETJjgHNfUAGPHJk44SocJPZ5sqYejjgp137nFhE7++OCD1Oe9xa0+/VS3zz7rbzypuJOwMcC6dYmPi+g481xK5DKhx9Pevbo9++xg48hReMffULTNm5f6vG2hL16sfywrVpQvJq9c+8vtjMCePYHDDtO+1LffTv9aFG1ffqnbQw8FVq2KTOscYEInv6Qbx21b6E8+qVvvhJ5ycn+Y/OUvyY+vXg2ccELiFO9zzwU2b2ZCjzM7ieiii/Q+i72XEgHsciF/DBuW+ry3y6VvX/9jSefNNzM/PmuWbu0iwHaIY6oPISb0eFi2zNl/6SXd9ukTTCwFYEInf7z/vm69l6orViQmv5qa8sVUKNtCz1QOIN+ZpRROf/5z8rmQrUqUCRM6+cvbcl20SEeTRIGd4m0TeqbFLOzNs717gXHjEodEUjTU1QG7diWeC/nMUC8mdCqNl15yZoUCTjL0ruhy/PHAypXOsU345V4sIJcuknPP1a2NrbY2/XNt6/2TT3RrL9cpOh56KPlcr15lD6MYWRO6iDwiIutF5NM0j4uI3Ccii0VkroicUPowKfTmzHHqtmzf7vSNewsZpZuZaRdbtrZvB2bPLi6eiRPTP26nc2di/5hbtdJWd7b1I5cudSYmRegynaC/b6lELKHnMsrlUQC/B/B4msfPBdC36WsYgD80balS3X134rE7ib/ySuJjtpKd97L2iSc0OfbrBxx0UP4xZGsh53IZ3bx5+sc6dgS2bHGOjQEed/2JeNcipXB74w3dXnyxdrHV1+txyBaBziZrQjfGTBGRXhmeciGAx40xBsB0EekgIocZY/gbTSpT90a6GaX51FEpRC5dLpkS+rXX6kSk7t2B3/xGhzK62fofFH5btgALFwKnngocd1zQ0RSlFH3o3QCsch3XNZ1LIiJjRWSmiMzc4K6ZQfFRV5d4vGhR5hEgNqF7W8w24RZ7QyrdB0axCb1tW70/0LIlhyxG3bJl+n84aFDQkRStrDdFjTETjDG1xpjazvmUJqXoSHVjyb2wrle6hG67YuzokUKl6yvPJQnb8eeFOPbYwr+XymfPHi1T0aaNdqNFXCkS+moAPVzH3ZvOEWVnE7eXTbi2NnlDAzBjRuahg6mka2XnktCLuToo5sOAyue114BNm4BLLw3los/5KsVPMAnAt5tGuwwHsI3955SzrVt1my552senTtU/vjlzSvO+NqFfcIFzrlvKnsLsLr00+VyIlykjODex580DTjkF6N076IhKIutvnYhMBDACQI2I1AH4FYDmAGCMGQ/gdQBjACwGsAvANX4FSzHSu3fiNOt0bIvcdp1s2aJLguU68iVdS9y+rvuDpHXr3F7TK1WdbLbQw+3ll3WYKRCbZA7kNsrlyiyPGwA/KllEVBnSJeQ9e5ySuoCusA44Ld4pU/Rr3Ljc3ifbQhvuhO5ekSaf+h1t2yaf++gjYPToWFzGx5J7clshw2JDiteFFIx0LdhXX9WE3ry5FvKyfeCF1kpJ932pErq7hT5kSO7vkS5p790bq2QRGw0NwM6dzvHBBwcXS4kxoVMw0g1TtH9oNhEXu/DF++8Dl13mHO/aBdx5p1MN0h1HqRcBZsGucFq+3LkZP2RI4V1tIcTrQSperuOw3SNO0t0E3bQpv9fMxlsh8Z57dPvhh7qtqnJmA9ruHSD/0TSpMKGH04IFugrVzTcn3hSPASZ0Kl6uyTdT98OaNcB776WvqZEv2w2ycKHGV1+v/e7eeuwiwHnn6WPuccj2gyVXqWYY3ntvfq9B5VFXpzeyYzgSiQmdipdra7ZrV2c/VQvdLshbCu6YjNHVh1JJd6UwY0Z+75euz93dV0vB27FDl5jr0SP7cyOICZ2Kl2k25zWuUazuFlE5h/U1Nqa/cZkuoee7klK6BMESF+Fii6ZFaBWifDChU/H+8Y/0j7nX7XT3oRd7uesedublvWIwJv+EfuGF+cVTVaWFurxSlR4YN85ZU5XKa8oU3XbpEmwcPmFCp/zNn5/YJeHtlwaAH/9YW+fuvuXaWqB9e90fNKi4CR2PPJL+MW88jY3JVwS2Be5N9LYfvZBp/3ZBDLd0pQ0WL87/9ak4+/c7xeO89fdjIn53Bch/dt1FOzok1WiOmprk9UJraoAbbnCOM1UzLIa3C+iZZ4Azz0w8Z0e/eBP3D39Y+AgX+/O0bAmMGaM12V98MfIlWWNh7VrgwQd1P9UHb0ywhU6FszXAc62cmW31olIwJnmBjWXLkocv2jU/vQm9efPCW282oTc0JC69t2ePs1YlS+0G47HHnP2jjw4uDp8xoVPh3n1Xt7YlbrsxBg9O/Xxv90Ypxnp72focXnZ5PMsOjyzlAsDuWa3uD6877tDJTBSM3bsT72WUe/3aMmJCp8LNnatbO5LDjhw4Ic2yst7k2alT6WNKt6CFbZF7+ZHQgcSEzlZ5MHbt0v/3VauyPzcm2IdOhaup0WT11lt63KcP8Mtf5l6QatQobdU//bQeV1fr5CO7nmMh3n47v+eXsniWO6Gz2mLwHngg8Xepffv8Ry9FDFvoVLhhw4Bbbkk8l0+CbNYMOOoo/RBo3VrrkefTDfPUU3o5PWmSM7Il1YibTErZQnf/7KV8Xcrfzp3JDYMbbojt+HOLLXQqnHdIXrpket55mWdMVlUBhx6q359P98QXX+gN0P37gcMO01E3LVvmt2xdqWdyjhiReTgmu1/8Zwzw/PPOcVUVcMklwcVTRkzoVLjJk3N7nh3emImI/iHme6PUfojY1nHnzvnVgyl2zVKvESOc/T590t+kBTLPYKXCff554uIpZ50FDBgQXDxlxN8mKp1Cl3ADNLE1NiYPL8zn+wFgyRLn3DHHpH6u+7wfI22sUaMyP/7FF/69dyWzyxaeeqpuhw4NLJRyY0Kn0imm31hECyelmqTkXkkonVQfBKefDhx+ePJ5d5eIn8uPecezjxvnjAwC0s8ipeLYq7bTTgN+9atYVlVMhwmd8jNrlj+vaxO62yGH6Na2dE89VRf0TeXVV5OXpauudvpO3YsYuEej2PfwQ6oJSi+/7OzzxmnpHTgAfPaZ/t83a1Zx/8aV89FFxXvqKf+6Cdwt7NGjtWVtSwx06wb8x384ibhrV+C557K/poizvJidqWnf67jj/G8hZ5tx6p49mqreDOXvzjv1/3XkyIpL5gATOuUjUzL/l38p7rXdxao6dtRRK7Z/u6oqsVWda199pj/oiy/OP8Z8ZUsoNqE/9ZT+/LkufE2pGeN8SJ98crCxBIRdLlS8448HvvKV0r2eXdvTndALYRNqly665JjlLRrml2wJ/d13dRw9Ky+Whp2xfN55FTt6qDJ/aspNY6Mu0pxutR+rmNEtqdj+7lGjdOZo27aJj7uTcya2G6d588TuFTsKIgzuuMPZ5xqkxfnkE03k/fsHHUlg2OVC6W3bpgvqrl0LXH99+d7XjkoYMCD1+OGDDtLx5tlWA7IJ0lum13vzNSz27auoERkltXIl8NFHwJFH5jYqKqbYQqf0cr2pVOqbT7kktWuvzf4cu1iF9woj02LVQcp3YWpSu3frgicHDpT+ajFimNCpMKNHO/v51k9JxT01O5fRHq1aOSNYUhk3zknc3vjCOvrBryGhcffpp85+BXe3AEzolIv6+uQaJO5W7scfF/8e/fqlfu1M0vWln3568fH4rVev5DHwdnk+yo87oZfrhndIMaFTdgcOAOPHJ55zr3J/2mnFv0chlQrTdc1kW5EmiPUkO3TQ6pTW17/udA+0a+c8h/LT0ACsW6f7rVuH9+qrTHgHhtJzt8rtH43lXpzCDjMshk3o+awzumVL6vOZFs644ALg2GNzf49iHXKI9o1fdhmwYoWeE9HkY0fe2Ju07lV1KDNjdF7E5s3673beeVy7FUzolEmupV5L1So6/3ztishVumF+3g+Y5s21H/3EE4EhQwoOryC2W0jE+dCy/67eDx4m9NxNngxMm+YcDx7MEUJglwtlkmslwlJN4qitza8PNNfKjHa0y8iR+cdULJtkDhzQIaBu3g+v3bvLElLkbd4MTJ/uHB95JJN5E/4rUHq5ttDXrg33SjBXXaX1sUvRNZQvO2KnoSF5XdOjjko8ZvXF3Eybpv+utsFx9tnBxhMibKFTeukS+tVXJx737Ol/LKnkenOzXbvg+lfdLfRLL8383EJrwVeatWu1QJvFm8n/hy10Si9Vgjn44OQa4qWs45KPKNTrOO00YPlyHdGS7QrBz8U24mLhQp0odtZZWpFTJJhRSyGVU0IXkdEA7gVQDeAhY8xvPI/3BPAYgA5Nz7nRGPN6aUOlsnvwwdyeF1TZ165dE1coAnQ4YJj07AncfLNzLKLrp6bCFnp6DQ26XOAzz+jxSSex3HAKWRO6iFQDuB/A2QDqAMwQkUnGmM9cT7sZwHPGmD+ISH8ArwPo5UO8FLTrrks+F1RL+bLLtLX2xBPOubDPFLz55sRRQddeCzz8sO7v2qVDMe1NXHLcdpuzf8klTOZp5PKXOBTAYmPMUmPMPgDPALjQ8xwDwM7Dbg9gTelCpFBx/yFdfrl+BTWZo1Ur4IgjnGM7QSfMqqsTPwDdE7RWrgTuvbf8MUVJdTUwcGDQUYRWLgm9G4BVruO6pnNu4wB8S0TqoK3zf031QiIyVkRmisjMDdkq5VHwTjop8+P9+iVO2Q9aptouFF3btzv7+Uw8q0Clula+EsCjxpjuAMYAeEJEkl7bGDPBGFNrjKnt3Llzid6afOMdPRDUaJZcRaGFTvm7+25nP4ihpxGSS0JfDcB1XYjuTefcrgXwHAAYYz4A0ApAZVfJiaLPPwfWr3eOvaMuwr6sV1jL4mZz6aW6qpL1pz/pwiKUPDZ/0KBAwoiKXBL6DAB9RaS3iLQAcAWASZ7nrAQwEgBEpB80obNPJWqefhp44AHn2D0O/ec/z170Kmh9+wYdQWGOPTZxKOiKFcmzSivVGs/tuFNOCSaOiMia0I0xBwD8GMCbABZAR7PMF5FbReSrTU/7fwC+JyJzAEwE8B1jcp1mSKFji1797W/OuTBPrbblAtyt3KjhWOrU5s939m+4Idy/hyGQ079O05jy1z3nfuna/wxAyK/HKWf33pu8An2u63gGwbYdolw6dfPm5HONjdGYPOUXY3Qi0eGHA+eey3rxOajg3xbKyH2Bdf754U6WcbgYTDV6o9Jru6xapaWFa2vTT8aiBEzopLxJ8aOPnP3a2vLGkq84tNBTLWxc6Ql9/nztYvEWMaO0mNAptbVrg44gd3YsfFRHuQCpu7T27i1/HGHy+ec6cYz3F3LGhE7AlCnJowmidIl71lnAL34R7YSeqmLgW2+VPYzQ2LZNb84ffnjQkUQKEzoBb78N/PGPief279ftd79b/njyVVWlS7pF2cCByXVoov4zFWrXLmDRIt1nd0temNAptb//XbesNV0eIomVIps3r7xZkfX1eqV4553A602D6vJZwYpYD52yaNs26Agqy6mnau2SlSs1wVWK3buB3/428VxQdfYjjC10yqySx0EHYeRI4KKLtIWebhHsuDFGyx14XXNN+WOJOLbQicJo3z6d/r9rV7z70uvrgblzE2sIWZXW5VQCbH5Vov37dZy5Mcnjz4cMCSYmSrR1q25ffTXQMHxljHazuEtMAMDo0cD3vhdMTBHHFnolmjoVePddHft8/PGJj7GmeLh89ln250SVu5pn167AmDH6+8ffwYKxhV6J7B/Stm3Jj+3eXd5YKLUzzgg6Av+57xE0awZ0785kXiQm9Epkb3ROnZq8MPH06eWPh5KdfnrQEfjPndBTlT6gvDGhVyKbxPftAz74IPEx96iWb36zfDFRenaSV9y4Sxtc6F2mmArBhF6JPvzQ2bcTiKyf/tTZ56SOcIhrQv/4Y91ecw1HtJQIE3olOuaY9I+5CyF17Oh/LJRdHMoDezU0AHPm6CpTrNdSMkzolShTESuuCBMep52m28ZGYM8e4LXXkouoRdW6dToj1jvKiorChF6JvIs/u0W5pnjcdOum261bgd/8BpgxA5gwIdCQSsZ+MPEqsKTYHKtE2S7hv/514JBDyhMLpWeT3apViefjMHt04ULdHnZYsHHEDBN6Jcq2cIK3jCsFw47J9hbp2rEj2gl9zx5g2TJg+HDWCioxJvRKlGqx3eOOA848s/yxUHr2BrX3A3jbtmhXIlywQG+KDhwYdCSxw4/HSpSqD/3001n7PGzs/YwFCxLPb99e/lhK6dNPgU6ddLo/lRQTeiVqaEgc93vmmewzD7Ndu3R79tnaRRHVhL55M1BXByxdChx7LG/A+4AJvRI1NADV1UCfPnocx3HOceHuWtm4EWjXLnUNnjAxBnjlFWD5cufczp3AffcBDz2kj9sRPFRSTOiVyPbJ2hIAdv1GCp9165x9Y/RGadhb6Fu26CzQRx91ir29+27icziZyBdM6JVg//7EVvi8edpismtYXnFFMHFRfvbs0YQe9ha6u1TBgw9qEa7Zs51zHTpwqr9PmNDjrqEBuP124LHHkh9r0wYYN44lS8PsvPOc/epqHaG0fXu4u8ncVRS3bgVWrEis6vmTn5Q9pEoRvWGLS5cCb70FXHVV5inspGz3iu3PDHMioGQnnKATi1q1AkaM0OXaDhzQsent2gUdXWretVBfeknj/8EP9PePY899E71/2RYtdNow+31z413xZvPmYOKgwlRXAxdfrKv5tG7tzCG4665g48rk2WcTj+vrdcx5hw6c6u+z6CX0bt00qX/5ZdCRRIN3OCKLb0Vb27bOvrclHAabNukwy3btgGHDnPODBwcXUwWJXkIX0YQ+a1bQkUSD+/J29erMhbko/NzdjGFcLtDWnbn66sRKipxEVBbRS+iAXsLt359c44KSufvMn3wS+Pzz4GKh4rmvsMK48MXKldpNdPDBugWAIUOCjamCRDOh2zGstmIbqZ07M99baNECeOMN3be1tila3DdCw5bQ9+0D5s/X/vIWLYAuXYBvfztxpA75KpoJ/eqrdetdD7PSPfccMHFi+kvxsLfuKLvqauBb39L9P/whXKOW5szRUVW2RS6is5E5qqVscvqXFpHRIrJIRBaLyI1pnvN1EflMROaLyNOlDdPD/oJs2uTr20TO1q263bPHOef+g3cndNahjq7mzZ1970zSoDQ2AtOna1959+7BxVHhsiZ0EakGcD+AcwH0B3CliPT3PKcvgJsAnGyMGQDg+tKH6nHWWbrdsMH3t4oMm7DT1Ttft067q9q0YenSKHOPbhk/Xrd/+xtwyy3ApEnBxLRkiTawhg9n0a0A5dJCHwpgsTFmqTFmH4BnAFzoec73ANxvjNkCAMaY9aUNMwXbCrj//sQWaSWzLfRMLbVNm/SGFf/oostb2KqhAZg2TfdnzwbWri1/TNOm6Rh5Lo4SqFwSejcA7jWw6prOuR0F4CgRmSoi00VkdKoXEpGxIjJTRGZuKLZl3aWLs7/e/8+PSLDTq90J3dufXl8fzB88lU6rVlqy4YIL9Ng72std5bAc9u7VsrhHH815DgEr1d2KZgD6AhgB4EoAfxSRDt4nGWMmGGNqjTG1nTt3Lu4dW7d2Js1wbHUim9BXrdIbpRRPti99/35dccoqd0mM2bM1Bve4cwpELgl9NYAeruPuTefc6gBMMsbsN8YsA/A5NMH765JLdPvOO76/Vei5+81tQn/44WBiofJo0UK3+/cnLis4bx7wxRfli2P6dL1iZo3zwOWS0GcA6CsivUWkBYArAHjvvLwMbZ1DRGqgXTBLSxdmGnbJNG81t0rk7mbJtgi0e5QERZddc9TbfblkCfDUU+W5cl20SMv5stszFLImdGPMAQA/BvAmgAUAnjPGzBeRW0Xkq01PexPAJhH5DMA7AH5ujPF/TGHr1k7p19tu8/3tIuOJJzJPMBo6tHyxkH/s7/6LL+oHup2ZafmdZPft05WJAGduCAUqpz50Y8zrxpijjDFHGGNubzr3S2PMpKZ9Y4z5qTGmvzFmoDHmGT+DTsC1MJV3ZMvEiYnHP/qRs3/GGf7HQ/7r1Em3NTXOuUMPdfbHjwfef78049ON0avgpUuBHTv03BtvaCGub38b6N27+PegokV/CldtrbNf6d0umbj7WL0tOYomEe273rgR+Mc/NOn+4AdAD9ctr7//vTQrHN1yi14FP/448MIL2m//8cf6WK9exb8+lUT0E/qAAcApp+j+p58GG0uQVnvvU3u0aAGcc45+cQx6fLi7VWyf+TnnJD6n2KGES5YkHi9frgttADrBj1P7QyMe/xP9+un2pZeCjSNIdlJRJv/0T/pF8ZGqy7FHDx2nXgoHDug9Ga+//EWvDk4+uTTvQyURj4Tu7kNc6v/gmlCaMSPoCCgIY8cC3/lO5uesWpX58Ux+/Wvd2iGSbkcfzau9kIlHQm/ZErjpJh2ON29e0NGUx/79ibNAbZGmVCvDnH56eWKi8mvZUvuwf/Yz4N//PfVzFiwo7LXdwx7/7d+SH3ffv6JQiEdCB/QXe+BA7UcP40oupXb77cAdd+iCAmvWOOePOirxeTfeyFEtlaBt2+RWtB3ZNHduYWPSb71Vt6NHaz95hw7aKu/eHRg5MvFGO4VCvAovDBqk05DvuEN/mYstLxAFjzySeNy9u87YW71aZ9K2ahVMXBQ89//9pEnAJ5/o/jnnZL+X4u7CO+YY7Vq57jp2sYRcfFroQOLU4/vv1xV8KknHjrqijR2mxhXWK1ubNs6+TeaAltr1XsW6ZxdPmQK89prun3qqMyObyTz04pXQq6sTh1D9z/+Ut6ZF0Ox44JEjdctJV5Wtqgr4z/9M/diuXc7+zJnAf/+31mTZvx94+209f8UV7K6LmHh1uQD6C7xtG/C73+nx/PlAX//rhJWFMdofmm6Iol1WbvDg1DdHqfJUV+sNzWnTdL7Ghx9qwl65Uh8/5BAtfQvoFe3Mmbo/Zox2tVCkxKuFDuhlYYcOTr2SjRsDDaekNm7UsfbpqktW8sQqSu+gg/SqrWVL56rtlVf0d8kY5yr2/feBN9/U/RNPDCZWKkr8Ero1ZoxubesjDrKN3hk0qCxhUIS51+Gtq9Pp/N57TQMGsL88ouKb0AGnGl25V3Dxw9SpTt9mOocfXp5YKLq8S8R5hzNefTVw2WXli4dKKn596G7DhgGTJ2uBomuvDTqawk2frj9HNtu3+x8LRVtNDXDhhcDmzdrFAgAXXaQrHrFVHnnxbqGffLJePq5fn33RhzDLlqhPOEG35V56jKJp8ODElvrxxzOZx0S8Ezqg/cp790Z7YWT3EDMgecr1iBHAVVcBQ4aULSSKuE6dgD59on3lSknin9Dt5IqottA/+ihxUgigkz3cC1a0bg0ccQTrnFPuWrbUhSnctdMp8ionoU+cqEtmRc3rryefa9ZMyxrYVWKKrXdNRLEQ/0xgR7oAwKOParnRqLOLPH/jG1qvmogIldBCFwG+9jXdX7NG10GMkrZtdXvSSTqcrLbWqarXvDlvhBLR/4l/Qgf0xqidxvzhh4GGkhdjgPp63R81SkfsnH9+sDERUWhVRkIHEuuEb94cXBzZHDigM/fmznVm9dlWOhFRBvHvQ7cGD9bx6NOnA/fdV7o1F0vtmWeAxYsTzw0bFkwsRBQpldNCF0lcDX3q1OBiycSbzAFO+iCinFROQge0PvSNN+r+5MlAQ4NWMDQm2Liyqaqs/yYiKkzlZYpWrZxRL7fdBvz+98Attzg3H8No2bKgIyCiCKi8hA5o7Qov79qcYcKbokSUg8pM6CJaJhRw1ksM88gXLjZARDmozIQO6LT5ceOA668POpLUzjzT2WcfOhHlgJnCLYhp9AsW6BqoXgMHOvsc5UJEOWBCB7TAPwA89VR53/f554FnnwXuuUdH2tg1QTt0ADp2dJ4XxaJiRFR2TOiArtbSvLmOJtm8uTwt9T17Ehd1XrJEE7yNB9Dp/gDQvr3/8RBR5DGhA9qlMXKk7t93H/DrX/v/nt6hiE8+6ezPm6fbk04CbropsWIkEVEaTOjW0KGJx95Vgkrt2WfTP+Ye1dKypb9xEFFs5JTQRWS0iCwSkcUicmOG510iIkZEatM9J7SqqpyhjEBxS9YVO/N0wIDivp+IKlLWhC4i1QDuB3AugP4ArhSR/ime1w7AdQAiVJ/Wo3dvZ2m3J54oLDEvWaIzT911143RlYeyfUhcd50uwME+cyIqQC4t9KEAFhtjlhpj9gF4BsCFKZ53G4A7AOwpYXzl5x5dcv/92vXS0ABs3Zr9e43RDwIgse76zp26NuiDDyY+f/jwxOMOHYCuXQuJmogop/K53QCsch3XAUio5yoiJwDoYYx5TUR+nu6FRGQsgLEA0LNnz/yjLYdmzYDLL9c+7o0bgTvvdB4bMya5r91t9+7kc7NnA5MmOccPPOCMopk+Hfj+951Ez/HmRFSEouuhi0gVgLsBfCfbc40xEwBMAIDa2trwljjs109HmHzwQeL511/XIl7DhwOtWyd/X2Ojs9+unSbsv/418Tnr1zv7AwcChx0GXHklx5oTUdFy6XJZDaCH67h70zmrHYBjAbwrIssBDAcwKZI3Rt1GjdLSAIcfnnh+yhTg7ru1C2bRIj3X2Ai89x6wZYvzvB07NJl37pz5PQDg6KMTZ4YSERUglxb6DAB9RaQ3NJFfAeAb9kFjzDYANfZYRN4F8DNjzMzShhqQq68GFi7Urpinn9ZzBw4Av/ud7p99tk5Keucd/fLasMF53uTJiY81q5wFo4jIf1kzijHmgIj8GMCbAKoBPGKMmS8itwKYaYyZlPkVIq6qCujfX5N4796avD//3Hncm6TT6dcv+bksukVEJSQmoNV6amtrzcyZEWzENzYCt96a3/eMGKFfxgDr1gHjx+v5m29mK52I8iIis4wxKbu0mU3yZVvV7dsDP/mJrnoEAJdeqkMcjzhC+9LdU/lHjNCtCHDooc756uqyhExElYEt9ELYfzMR4JNPtNZKnz7Jz3v/faBLF73p6TZvnn6P94YrEVEWbKGXmnu8+KBB6Z936qmpz3NECxH5gHfliIhiggmdiCgmmNCJiGKCCZ2IKCaY0ImIYoIJnYgoJpjQiYhiggmdiCgmApspKiIbAKwo8NtrAGwsYThB4s8STnH5WeLycwD8WazDjTEp63IHltCLISIz0019jRr+LOEUl58lLj8HwJ8lF+xyISKKCSZ0IqKYiGpCnxB0ACXEnyWc4vKzxOXnAPizZBXJPnQiIkoW1RY6ERF5MKETEcVEpBK6iDwiIutF5NOgYymWiPQQkXdE5DMRmS8i1wUdUyFEpJWIfCQic5p+jluCjqlYIlItIh+LyKtBx1IMEVkuIvNE5BMRiejyYEpEOojI8yKyUEQWiMhJQceULxE5uun/wn5tF5HrS/oeUepDF5HTANQDeNwYc2zQ8RRDRA4DcJgxZraItAMwC8DXjDGfBRxaXkREALQxxtSLSHMA/wBwnTFmesChFUxEfgqgFsDBxpjzg46nUCKyHECtMSbyk3FE5DEA7xtjHhKRFgBaG2O2BhxWwUSkGsBqAMOMMYVOsEwSqRa6MWYKgM1Bx1EKxpi1xpjZTfs7ACwA0C3YqPJnVH3TYfOmr+i0EjxEpDuA8wA8FHQspESkPYDTADwMAMaYfVFO5k1GAlhSymQORCyhx5WI9AIwGMCHAYdSkKYuik8ArAcw2RgTyZ+jye8A/AJAY8BxlIIB8DcRmSUiY4MOpgi9AWwA8KemrrCHRKRN0EEV6QoAE0v9okzoARORtgBeAHC9MWZ70PEUwhjTYIwZBKA7gKEiEsnuMBE5H8B6Y8ysoGMpkVOMMScAOBfAj5q6LKOoGYATAPzBGDMYwE4ANwYbUuGauoy+CuDPpX5tJvQANfU5vwDgKWPMi0HHU6ymy+B3AIwOOJRCnQzgq019z88AOFNEngw2pMIZY1Y3bdcDeAnA0GAjKlgdgDrXld/z0AQfVecCmG2MWVfqF2ZCD0jTzcSHASwwxtwddDyFEpHOItKhaf8gAGcDWBhoUAUyxtxkjOlujOkFvSR+2xjzrYDDKoiItGm62Y6m7olzAERydJgx5ksAq0Tk6KZTIwFEavCAx5XwobsF0EuZyBCRiQBGAKgRkToAvzLGPBxsVAU7GcBVAOY19T8DwL8bY14PLqSCHAbgsaa79lUAnjPGRHq4X0x8BcBL2m5AMwBPG2P+GmxIRflXAE81dVcsBXBNwPEUpOnD9WwA3/fl9aM0bJGIiNJjlwsRUUwwoRMRxQQTOhFRTDChExHFBBM6EVFMMKETEcUEEzoRUUz8fyocXx+6ABl8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "out = []\n",
    "rs = np.linspace(1,7,10000)\n",
    "for r in tqdm(rs):\n",
    "    out.append(goal_function(r,2,2))\n",
    "plt.plot(rs,out, c='r', alpha=.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "966512c4",
   "metadata": {},
   "source": [
    "## Free param mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 479,
   "id": "8fdfb1ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "rho_c_for_mesh = np.linspace(3.4,3.8,20)\n",
    "nu_for_mesh = np.linspace(1.95,2.05,1)\n",
    "zeta_for_mesh = np.linspace(1.95,2.05,1)\n",
    "\n",
    "rho_c_mesh, nu_mesh, zeta_mesh = np.meshgrid(rho_c_for_mesh,nu_for_mesh,zeta_for_mesh,indexing='ij')\n",
    "quality_mesh = np.zeros(np.shape(rho_c_mesh))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "624ace53",
   "metadata": {},
   "source": [
    "## calculate quality for each mesh point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 480,
   "id": "a8b50f74",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "20it [00:00, 72.46it/s]\n"
     ]
    }
   ],
   "source": [
    "real_data = []\n",
    "\n",
    "for index3, rho_3, nu_3, z_3 in tqdm(zip(range(len(rho_c_mesh)),rho_c_mesh , nu_mesh, zeta_mesh)):\n",
    "    for index2, rho_2, nu_2, z_2 in zip(range(len(rho_3)),rho_3, nu_3, z_3):\n",
    "        for index1, rho_1, nu_1, z_1 in zip(range(len(rho_2)),rho_2, nu_2, z_2):\n",
    "            q = goal_function(rho_1, nu_1, z_1)\n",
    "            quality_mesh[index3,index2,index1] = q\n",
    "            real_data.append([rho_1, nu_1, z_1, q])\n",
    "            #print(quality_mesh)\n",
    "            \n",
    "real_data = np.array(real_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 481,
   "id": "5624afa7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 481,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(real_data[:,3] != quality_mesh.flatten())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "949e7a81",
   "metadata": {},
   "source": [
    "## Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 482,
   "id": "edeff3c4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.plotly.v1+json": {
       "config": {
        "plotlyServerURL": "https://plot.ly"
       },
       "data": [
        {
         "marker": {
          "color": [
           5.375960851644816,
           8.406103005113401,
           6.358585851482232,
           6.321442175804631,
           5.523292438037945,
           7.794604824431682,
           7.6059427972017035,
           10.088257820332334,
           10.487998550743288,
           12.722867658880224,
           14.760268210596534,
           11.266801774141971,
           13.21156548105472,
           9.71264732679555,
           9.124532163477488,
           4.718359108512384,
           6.230564344463909,
           7.54721901548738,
           11.027831345405865,
           6.717108443277575
          ],
          "colorscale": [
           [
            0,
            "#000004"
           ],
           [
            0.1111111111111111,
            "#1b0c41"
           ],
           [
            0.2222222222222222,
            "#4a0c6b"
           ],
           [
            0.3333333333333333,
            "#781c6d"
           ],
           [
            0.4444444444444444,
            "#a52c60"
           ],
           [
            0.5555555555555556,
            "#cf4446"
           ],
           [
            0.6666666666666666,
            "#ed6925"
           ],
           [
            0.7777777777777778,
            "#fb9b06"
           ],
           [
            0.8888888888888888,
            "#f7d13d"
           ],
           [
            1,
            "#fcffa4"
           ]
          ],
          "opacity": 0.8,
          "size": [
           5.375960851644816,
           8.406103005113401,
           6.358585851482232,
           6.321442175804631,
           5.523292438037945,
           7.794604824431682,
           7.6059427972017035,
           10.088257820332334,
           10.487998550743288,
           12.722867658880224,
           14.760268210596534,
           11.266801774141971,
           13.21156548105472,
           9.71264732679555,
           9.124532163477488,
           4.718359108512384,
           6.230564344463909,
           7.54721901548738,
           11.027831345405865,
           6.717108443277575
          ]
         },
         "mode": "markers",
         "type": "scatter3d",
         "x": [
          3.4,
          3.4210526315789473,
          3.442105263157895,
          3.463157894736842,
          3.484210526315789,
          3.5052631578947366,
          3.526315789473684,
          3.5473684210526315,
          3.568421052631579,
          3.5894736842105264,
          3.6105263157894734,
          3.631578947368421,
          3.6526315789473682,
          3.6736842105263157,
          3.694736842105263,
          3.7157894736842105,
          3.7368421052631575,
          3.757894736842105,
          3.7789473684210524,
          3.8
         ],
         "y": [
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95
         ],
         "z": [
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95,
          1.95
         ]
        }
       ],
       "layout": {
        "template": {
         "data": {
          "bar": [
           {
            "error_x": {
             "color": "#2a3f5f"
            },
            "error_y": {
             "color": "#2a3f5f"
            },
            "marker": {
             "line": {
              "color": "#E5ECF6",
              "width": 0.5
             },
             "pattern": {
              "fillmode": "overlay",
              "size": 10,
              "solidity": 0.2
             }
            },
            "type": "bar"
           }
          ],
          "barpolar": [
           {
            "marker": {
             "line": {
              "color": "#E5ECF6",
              "width": 0.5
             },
             "pattern": {
              "fillmode": "overlay",
              "size": 10,
              "solidity": 0.2
             }
            },
            "type": "barpolar"
           }
          ],
          "carpet": [
           {
            "aaxis": {
             "endlinecolor": "#2a3f5f",
             "gridcolor": "white",
             "linecolor": "white",
             "minorgridcolor": "white",
             "startlinecolor": "#2a3f5f"
            },
            "baxis": {
             "endlinecolor": "#2a3f5f",
             "gridcolor": "white",
             "linecolor": "white",
             "minorgridcolor": "white",
             "startlinecolor": "#2a3f5f"
            },
            "type": "carpet"
           }
          ],
          "choropleth": [
           {
            "colorbar": {
             "outlinewidth": 0,
             "ticks": ""
            },
            "type": "choropleth"
           }
          ],
          "contour": [
           {
            "colorbar": {
             "outlinewidth": 0,
             "ticks": ""
            },
            "colorscale": [
             [
              0,
              "#0d0887"
             ],
             [
              0.1111111111111111,
              "#46039f"
             ],
             [
              0.2222222222222222,
              "#7201a8"
             ],
             [
              0.3333333333333333,
              "#9c179e"
             ],
             [
              0.4444444444444444,
              "#bd3786"
             ],
             [
              0.5555555555555556,
              "#d8576b"
             ],
             [
              0.6666666666666666,
              "#ed7953"
             ],
             [
              0.7777777777777778,
              "#fb9f3a"
             ],
             [
              0.8888888888888888,
              "#fdca26"
             ],
             [
              1,
              "#f0f921"
             ]
            ],
            "type": "contour"
           }
          ],
          "contourcarpet": [
           {
            "colorbar": {
             "outlinewidth": 0,
             "ticks": ""
            },
            "type": "contourcarpet"
           }
          ],
          "heatmap": [
           {
            "colorbar": {
             "outlinewidth": 0,
             "ticks": ""
            },
            "colorscale": [
             [
              0,
              "#0d0887"
             ],
             [
              0.1111111111111111,
              "#46039f"
             ],
             [
              0.2222222222222222,
              "#7201a8"
             ],
             [
              0.3333333333333333,
              "#9c179e"
             ],
             [
              0.4444444444444444,
              "#bd3786"
             ],
             [
              0.5555555555555556,
              "#d8576b"
             ],
             [
              0.6666666666666666,
              "#ed7953"
             ],
             [
              0.7777777777777778,
              "#fb9f3a"
             ],
             [
              0.8888888888888888,
              "#fdca26"
             ],
             [
              1,
              "#f0f921"
             ]
            ],
            "type": "heatmap"
           }
          ],
          "heatmapgl": [
           {
            "colorbar": {
             "outlinewidth": 0,
             "ticks": ""
            },
            "colorscale": [
             [
              0,
              "#0d0887"
             ],
             [
              0.1111111111111111,
              "#46039f"
             ],
             [
              0.2222222222222222,
              "#7201a8"
             ],
             [
              0.3333333333333333,
              "#9c179e"
             ],
             [
              0.4444444444444444,
              "#bd3786"
             ],
             [
              0.5555555555555556,
              "#d8576b"
             ],
             [
              0.6666666666666666,
              "#ed7953"
             ],
             [
              0.7777777777777778,
              "#fb9f3a"
             ],
             [
              0.8888888888888888,
              "#fdca26"
             ],
             [
              1,
              "#f0f921"
             ]
            ],
            "type": "heatmapgl"
           }
          ],
          "histogram": [
           {
            "marker": {
             "pattern": {
              "fillmode": "overlay",
              "size": 10,
              "solidity": 0.2
             }
            },
            "type": "histogram"
           }
          ],
          "histogram2d": [
           {
            "colorbar": {
             "outlinewidth": 0,
             "ticks": ""
            },
            "colorscale": [
             [
              0,
              "#0d0887"
             ],
             [
              0.1111111111111111,
              "#46039f"
             ],
             [
              0.2222222222222222,
              "#7201a8"
             ],
             [
              0.3333333333333333,
              "#9c179e"
             ],
             [
              0.4444444444444444,
              "#bd3786"
             ],
             [
              0.5555555555555556,
              "#d8576b"
             ],
             [
              0.6666666666666666,
              "#ed7953"
             ],
             [
              0.7777777777777778,
              "#fb9f3a"
             ],
             [
              0.8888888888888888,
              "#fdca26"
             ],
             [
              1,
              "#f0f921"
             ]
            ],
            "type": "histogram2d"
           }
          ],
          "histogram2dcontour": [
           {
            "colorbar": {
             "outlinewidth": 0,
             "ticks": ""
            },
            "colorscale": [
             [
              0,
              "#0d0887"
             ],
             [
              0.1111111111111111,
              "#46039f"
             ],
             [
              0.2222222222222222,
              "#7201a8"
             ],
             [
              0.3333333333333333,
              "#9c179e"
             ],
             [
              0.4444444444444444,
              "#bd3786"
             ],
             [
              0.5555555555555556,
              "#d8576b"
             ],
             [
              0.6666666666666666,
              "#ed7953"
             ],
             [
              0.7777777777777778,
              "#fb9f3a"
             ],
             [
              0.8888888888888888,
              "#fdca26"
             ],
             [
              1,
              "#f0f921"
             ]
            ],
            "type": "histogram2dcontour"
           }
          ],
          "mesh3d": [
           {
            "colorbar": {
             "outlinewidth": 0,
             "ticks": ""
            },
            "type": "mesh3d"
           }
          ],
          "parcoords": [
           {
            "line": {
             "colorbar": {
              "outlinewidth": 0,
              "ticks": ""
             }
            },
            "type": "parcoords"
           }
          ],
          "pie": [
           {
            "automargin": true,
            "type": "pie"
           }
          ],
          "scatter": [
           {
            "marker": {
             "colorbar": {
              "outlinewidth": 0,
              "ticks": ""
             }
            },
            "type": "scatter"
           }
          ],
          "scatter3d": [
           {
            "line": {
             "colorbar": {
              "outlinewidth": 0,
              "ticks": ""
             }
            },
            "marker": {
             "colorbar": {
              "outlinewidth": 0,
              "ticks": ""
             }
            },
            "type": "scatter3d"
           }
          ],
          "scattercarpet": [
           {
            "marker": {
             "colorbar": {
              "outlinewidth": 0,
              "ticks": ""
             }
            },
            "type": "scattercarpet"
           }
          ],
          "scattergeo": [
           {
            "marker": {
             "colorbar": {
              "outlinewidth": 0,
              "ticks": ""
             }
            },
            "type": "scattergeo"
           }
          ],
          "scattergl": [
           {
            "marker": {
             "colorbar": {
              "outlinewidth": 0,
              "ticks": ""
             }
            },
            "type": "scattergl"
           }
          ],
          "scattermapbox": [
           {
            "marker": {
             "colorbar": {
              "outlinewidth": 0,
              "ticks": ""
             }
            },
            "type": "scattermapbox"
           }
          ],
          "scatterpolar": [
           {
            "marker": {
             "colorbar": {
              "outlinewidth": 0,
              "ticks": ""
             }
            },
            "type": "scatterpolar"
           }
          ],
          "scatterpolargl": [
           {
            "marker": {
             "colorbar": {
              "outlinewidth": 0,
              "ticks": ""
             }
            },
            "type": "scatterpolargl"
           }
          ],
          "scatterternary": [
           {
            "marker": {
             "colorbar": {
              "outlinewidth": 0,
              "ticks": ""
             }
            },
            "type": "scatterternary"
           }
          ],
          "surface": [
           {
            "colorbar": {
             "outlinewidth": 0,
             "ticks": ""
            },
            "colorscale": [
             [
              0,
              "#0d0887"
             ],
             [
              0.1111111111111111,
              "#46039f"
             ],
             [
              0.2222222222222222,
              "#7201a8"
             ],
             [
              0.3333333333333333,
              "#9c179e"
             ],
             [
              0.4444444444444444,
              "#bd3786"
             ],
             [
              0.5555555555555556,
              "#d8576b"
             ],
             [
              0.6666666666666666,
              "#ed7953"
             ],
             [
              0.7777777777777778,
              "#fb9f3a"
             ],
             [
              0.8888888888888888,
              "#fdca26"
             ],
             [
              1,
              "#f0f921"
             ]
            ],
            "type": "surface"
           }
          ],
          "table": [
           {
            "cells": {
             "fill": {
              "color": "#EBF0F8"
             },
             "line": {
              "color": "white"
             }
            },
            "header": {
             "fill": {
              "color": "#C8D4E3"
             },
             "line": {
              "color": "white"
             }
            },
            "type": "table"
           }
          ]
         },
         "layout": {
          "annotationdefaults": {
           "arrowcolor": "#2a3f5f",
           "arrowhead": 0,
           "arrowwidth": 1
          },
          "autotypenumbers": "strict",
          "coloraxis": {
           "colorbar": {
            "outlinewidth": 0,
            "ticks": ""
           }
          },
          "colorscale": {
           "diverging": [
            [
             0,
             "#8e0152"
            ],
            [
             0.1,
             "#c51b7d"
            ],
            [
             0.2,
             "#de77ae"
            ],
            [
             0.3,
             "#f1b6da"
            ],
            [
             0.4,
             "#fde0ef"
            ],
            [
             0.5,
             "#f7f7f7"
            ],
            [
             0.6,
             "#e6f5d0"
            ],
            [
             0.7,
             "#b8e186"
            ],
            [
             0.8,
             "#7fbc41"
            ],
            [
             0.9,
             "#4d9221"
            ],
            [
             1,
             "#276419"
            ]
           ],
           "sequential": [
            [
             0,
             "#0d0887"
            ],
            [
             0.1111111111111111,
             "#46039f"
            ],
            [
             0.2222222222222222,
             "#7201a8"
            ],
            [
             0.3333333333333333,
             "#9c179e"
            ],
            [
             0.4444444444444444,
             "#bd3786"
            ],
            [
             0.5555555555555556,
             "#d8576b"
            ],
            [
             0.6666666666666666,
             "#ed7953"
            ],
            [
             0.7777777777777778,
             "#fb9f3a"
            ],
            [
             0.8888888888888888,
             "#fdca26"
            ],
            [
             1,
             "#f0f921"
            ]
           ],
           "sequentialminus": [
            [
             0,
             "#0d0887"
            ],
            [
             0.1111111111111111,
             "#46039f"
            ],
            [
             0.2222222222222222,
             "#7201a8"
            ],
            [
             0.3333333333333333,
             "#9c179e"
            ],
            [
             0.4444444444444444,
             "#bd3786"
            ],
            [
             0.5555555555555556,
             "#d8576b"
            ],
            [
             0.6666666666666666,
             "#ed7953"
            ],
            [
             0.7777777777777778,
             "#fb9f3a"
            ],
            [
             0.8888888888888888,
             "#fdca26"
            ],
            [
             1,
             "#f0f921"
            ]
           ]
          },
          "colorway": [
           "#636efa",
           "#EF553B",
           "#00cc96",
           "#ab63fa",
           "#FFA15A",
           "#19d3f3",
           "#FF6692",
           "#B6E880",
           "#FF97FF",
           "#FECB52"
          ],
          "font": {
           "color": "#2a3f5f"
          },
          "geo": {
           "bgcolor": "white",
           "lakecolor": "white",
           "landcolor": "#E5ECF6",
           "showlakes": true,
           "showland": true,
           "subunitcolor": "white"
          },
          "hoverlabel": {
           "align": "left"
          },
          "hovermode": "closest",
          "mapbox": {
           "style": "light"
          },
          "paper_bgcolor": "white",
          "plot_bgcolor": "#E5ECF6",
          "polar": {
           "angularaxis": {
            "gridcolor": "white",
            "linecolor": "white",
            "ticks": ""
           },
           "bgcolor": "#E5ECF6",
           "radialaxis": {
            "gridcolor": "white",
            "linecolor": "white",
            "ticks": ""
           }
          },
          "scene": {
           "xaxis": {
            "backgroundcolor": "#E5ECF6",
            "gridcolor": "white",
            "gridwidth": 2,
            "linecolor": "white",
            "showbackground": true,
            "ticks": "",
            "zerolinecolor": "white"
           },
           "yaxis": {
            "backgroundcolor": "#E5ECF6",
            "gridcolor": "white",
            "gridwidth": 2,
            "linecolor": "white",
            "showbackground": true,
            "ticks": "",
            "zerolinecolor": "white"
           },
           "zaxis": {
            "backgroundcolor": "#E5ECF6",
            "gridcolor": "white",
            "gridwidth": 2,
            "linecolor": "white",
            "showbackground": true,
            "ticks": "",
            "zerolinecolor": "white"
           }
          },
          "shapedefaults": {
           "line": {
            "color": "#2a3f5f"
           }
          },
          "ternary": {
           "aaxis": {
            "gridcolor": "white",
            "linecolor": "white",
            "ticks": ""
           },
           "baxis": {
            "gridcolor": "white",
            "linecolor": "white",
            "ticks": ""
           },
           "bgcolor": "#E5ECF6",
           "caxis": {
            "gridcolor": "white",
            "linecolor": "white",
            "ticks": ""
           }
          },
          "title": {
           "x": 0.05
          },
          "xaxis": {
           "automargin": true,
           "gridcolor": "white",
           "linecolor": "white",
           "ticks": "",
           "title": {
            "standoff": 15
           },
           "zerolinecolor": "white",
           "zerolinewidth": 2
          },
          "yaxis": {
           "automargin": true,
           "gridcolor": "white",
           "linecolor": "white",
           "ticks": "",
           "title": {
            "standoff": 15
           },
           "zerolinecolor": "white",
           "zerolinewidth": 2
          }
         }
        }
       }
      },
      "text/html": [
       "<div>                            <div id=\"5b0d60c3-9682-4d91-a2dd-c0e82e7fcd49\" class=\"plotly-graph-div\" style=\"height:525px; width:100%;\"></div>            <script type=\"text/javascript\">                require([\"plotly\"], function(Plotly) {                    window.PLOTLYENV=window.PLOTLYENV || {};                                    if (document.getElementById(\"5b0d60c3-9682-4d91-a2dd-c0e82e7fcd49\")) {                    Plotly.newPlot(                        \"5b0d60c3-9682-4d91-a2dd-c0e82e7fcd49\",                        [{\"marker\":{\"color\":[5.375960851644816,8.406103005113401,6.358585851482232,6.321442175804631,5.523292438037945,7.794604824431682,7.6059427972017035,10.088257820332334,10.487998550743288,12.722867658880224,14.760268210596534,11.266801774141971,13.21156548105472,9.71264732679555,9.124532163477488,4.718359108512384,6.230564344463909,7.54721901548738,11.027831345405865,6.717108443277575],\"colorscale\":[[0.0,\"#000004\"],[0.1111111111111111,\"#1b0c41\"],[0.2222222222222222,\"#4a0c6b\"],[0.3333333333333333,\"#781c6d\"],[0.4444444444444444,\"#a52c60\"],[0.5555555555555556,\"#cf4446\"],[0.6666666666666666,\"#ed6925\"],[0.7777777777777778,\"#fb9b06\"],[0.8888888888888888,\"#f7d13d\"],[1.0,\"#fcffa4\"]],\"opacity\":0.8,\"size\":[5.375960851644816,8.406103005113401,6.358585851482232,6.321442175804631,5.523292438037945,7.794604824431682,7.6059427972017035,10.088257820332334,10.487998550743288,12.722867658880224,14.760268210596534,11.266801774141971,13.21156548105472,9.71264732679555,9.124532163477488,4.718359108512384,6.230564344463909,7.54721901548738,11.027831345405865,6.717108443277575]},\"mode\":\"markers\",\"type\":\"scatter3d\",\"x\":[3.4,3.4210526315789473,3.442105263157895,3.463157894736842,3.484210526315789,3.5052631578947366,3.526315789473684,3.5473684210526315,3.568421052631579,3.5894736842105264,3.6105263157894734,3.631578947368421,3.6526315789473682,3.6736842105263157,3.694736842105263,3.7157894736842105,3.7368421052631575,3.757894736842105,3.7789473684210524,3.8],\"y\":[1.95,1.95,1.95,1.95,1.95,1.95,1.95,1.95,1.95,1.95,1.95,1.95,1.95,1.95,1.95,1.95,1.95,1.95,1.95,1.95],\"z\":[1.95,1.95,1.95,1.95,1.95,1.95,1.95,1.95,1.95,1.95,1.95,1.95,1.95,1.95,1.95,1.95,1.95,1.95,1.95,1.95]}],                        {\"template\":{\"data\":{\"bar\":[{\"error_x\":{\"color\":\"#2a3f5f\"},\"error_y\":{\"color\":\"#2a3f5f\"},\"marker\":{\"line\":{\"color\":\"#E5ECF6\",\"width\":0.5},\"pattern\":{\"fillmode\":\"overlay\",\"size\":10,\"solidity\":0.2}},\"type\":\"bar\"}],\"barpolar\":[{\"marker\":{\"line\":{\"color\":\"#E5ECF6\",\"width\":0.5},\"pattern\":{\"fillmode\":\"overlay\",\"size\":10,\"solidity\":0.2}},\"type\":\"barpolar\"}],\"carpet\":[{\"aaxis\":{\"endlinecolor\":\"#2a3f5f\",\"gridcolor\":\"white\",\"linecolor\":\"white\",\"minorgridcolor\":\"white\",\"startlinecolor\":\"#2a3f5f\"},\"baxis\":{\"endlinecolor\":\"#2a3f5f\",\"gridcolor\":\"white\",\"linecolor\":\"white\",\"minorgridcolor\":\"white\",\"startlinecolor\":\"#2a3f5f\"},\"type\":\"carpet\"}],\"choropleth\":[{\"colorbar\":{\"outlinewidth\":0,\"ticks\":\"\"},\"type\":\"choropleth\"}],\"contour\":[{\"colorbar\":{\"outlinewidth\":0,\"ticks\":\"\"},\"colorscale\":[[0.0,\"#0d0887\"],[0.1111111111111111,\"#46039f\"],[0.2222222222222222,\"#7201a8\"],[0.3333333333333333,\"#9c179e\"],[0.4444444444444444,\"#bd3786\"],[0.5555555555555556,\"#d8576b\"],[0.6666666666666666,\"#ed7953\"],[0.7777777777777778,\"#fb9f3a\"],[0.8888888888888888,\"#fdca26\"],[1.0,\"#f0f921\"]],\"type\":\"contour\"}],\"contourcarpet\":[{\"colorbar\":{\"outlinewidth\":0,\"ticks\":\"\"},\"type\":\"contourcarpet\"}],\"heatmap\":[{\"colorbar\":{\"outlinewidth\":0,\"ticks\":\"\"},\"colorscale\":[[0.0,\"#0d0887\"],[0.1111111111111111,\"#46039f\"],[0.2222222222222222,\"#7201a8\"],[0.3333333333333333,\"#9c179e\"],[0.4444444444444444,\"#bd3786\"],[0.5555555555555556,\"#d8576b\"],[0.6666666666666666,\"#ed7953\"],[0.7777777777777778,\"#fb9f3a\"],[0.8888888888888888,\"#fdca26\"],[1.0,\"#f0f921\"]],\"type\":\"heatmap\"}],\"heatmapgl\":[{\"colorbar\":{\"outlinewidth\":0,\"ticks\":\"\"},\"colorscale\":[[0.0,\"#0d0887\"],[0.1111111111111111,\"#46039f\"],[0.2222222222222222,\"#7201a8\"],[0.3333333333333333,\"#9c179e\"],[0.4444444444444444,\"#bd3786\"],[0.5555555555555556,\"#d8576b\"],[0.6666666666666666,\"#ed7953\"],[0.7777777777777778,\"#fb9f3a\"],[0.8888888888888888,\"#fdca26\"],[1.0,\"#f0f921\"]],\"type\":\"heatmapgl\"}],\"histogram\":[{\"marker\":{\"pattern\":{\"fillmode\":\"overlay\",\"size\":10,\"solidity\":0.2}},\"type\":\"histogram\"}],\"histogram2d\":[{\"colorbar\":{\"outlinewidth\":0,\"ticks\":\"\"},\"colorscale\":[[0.0,\"#0d0887\"],[0.1111111111111111,\"#46039f\"],[0.2222222222222222,\"#7201a8\"],[0.3333333333333333,\"#9c179e\"],[0.4444444444444444,\"#bd3786\"],[0.5555555555555556,\"#d8576b\"],[0.6666666666666666,\"#ed7953\"],[0.7777777777777778,\"#fb9f3a\"],[0.8888888888888888,\"#fdca26\"],[1.0,\"#f0f921\"]],\"type\":\"histogram2d\"}],\"histogram2dcontour\":[{\"colorbar\":{\"outlinewidth\":0,\"ticks\":\"\"},\"colorscale\":[[0.0,\"#0d0887\"],[0.1111111111111111,\"#46039f\"],[0.2222222222222222,\"#7201a8\"],[0.3333333333333333,\"#9c179e\"],[0.4444444444444444,\"#bd3786\"],[0.5555555555555556,\"#d8576b\"],[0.6666666666666666,\"#ed7953\"],[0.7777777777777778,\"#fb9f3a\"],[0.8888888888888888,\"#fdca26\"],[1.0,\"#f0f921\"]],\"type\":\"histogram2dcontour\"}],\"mesh3d\":[{\"colorbar\":{\"outlinewidth\":0,\"ticks\":\"\"},\"type\":\"mesh3d\"}],\"parcoords\":[{\"line\":{\"colorbar\":{\"outlinewidth\":0,\"ticks\":\"\"}},\"type\":\"parcoords\"}],\"pie\":[{\"automargin\":true,\"type\":\"pie\"}],\"scatter\":[{\"marker\":{\"colorbar\":{\"outlinewidth\":0,\"ticks\":\"\"}},\"type\":\"scatter\"}],\"scatter3d\":[{\"line\":{\"colorbar\":{\"outlinewidth\":0,\"ticks\":\"\"}},\"marker\":{\"colorbar\":{\"outlinewidth\":0,\"ticks\":\"\"}},\"type\":\"scatter3d\"}],\"scattercarpet\":[{\"marker\":{\"colorbar\":{\"outlinewidth\":0,\"ticks\":\"\"}},\"type\":\"scattercarpet\"}],\"scattergeo\":[{\"marker\":{\"colorbar\":{\"outlinewidth\":0,\"ticks\":\"\"}},\"type\":\"scattergeo\"}],\"scattergl\":[{\"marker\":{\"colorbar\":{\"outlinewidth\":0,\"ticks\":\"\"}},\"type\":\"scattergl\"}],\"scattermapbox\":[{\"marker\":{\"colorbar\":{\"outlinewidth\":0,\"ticks\":\"\"}},\"type\":\"scattermapbox\"}],\"scatterpolar\":[{\"marker\":{\"colorbar\":{\"outlinewidth\":0,\"ticks\":\"\"}},\"type\":\"scatterpolar\"}],\"scatterpolargl\":[{\"marker\":{\"colorbar\":{\"outlinewidth\":0,\"ticks\":\"\"}},\"type\":\"scatterpolargl\"}],\"scatterternary\":[{\"marker\":{\"colorbar\":{\"outlinewidth\":0,\"ticks\":\"\"}},\"type\":\"scatterternary\"}],\"surface\":[{\"colorbar\":{\"outlinewidth\":0,\"ticks\":\"\"},\"colorscale\":[[0.0,\"#0d0887\"],[0.1111111111111111,\"#46039f\"],[0.2222222222222222,\"#7201a8\"],[0.3333333333333333,\"#9c179e\"],[0.4444444444444444,\"#bd3786\"],[0.5555555555555556,\"#d8576b\"],[0.6666666666666666,\"#ed7953\"],[0.7777777777777778,\"#fb9f3a\"],[0.8888888888888888,\"#fdca26\"],[1.0,\"#f0f921\"]],\"type\":\"surface\"}],\"table\":[{\"cells\":{\"fill\":{\"color\":\"#EBF0F8\"},\"line\":{\"color\":\"white\"}},\"header\":{\"fill\":{\"color\":\"#C8D4E3\"},\"line\":{\"color\":\"white\"}},\"type\":\"table\"}]},\"layout\":{\"annotationdefaults\":{\"arrowcolor\":\"#2a3f5f\",\"arrowhead\":0,\"arrowwidth\":1},\"autotypenumbers\":\"strict\",\"coloraxis\":{\"colorbar\":{\"outlinewidth\":0,\"ticks\":\"\"}},\"colorscale\":{\"diverging\":[[0,\"#8e0152\"],[0.1,\"#c51b7d\"],[0.2,\"#de77ae\"],[0.3,\"#f1b6da\"],[0.4,\"#fde0ef\"],[0.5,\"#f7f7f7\"],[0.6,\"#e6f5d0\"],[0.7,\"#b8e186\"],[0.8,\"#7fbc41\"],[0.9,\"#4d9221\"],[1,\"#276419\"]],\"sequential\":[[0.0,\"#0d0887\"],[0.1111111111111111,\"#46039f\"],[0.2222222222222222,\"#7201a8\"],[0.3333333333333333,\"#9c179e\"],[0.4444444444444444,\"#bd3786\"],[0.5555555555555556,\"#d8576b\"],[0.6666666666666666,\"#ed7953\"],[0.7777777777777778,\"#fb9f3a\"],[0.8888888888888888,\"#fdca26\"],[1.0,\"#f0f921\"]],\"sequentialminus\":[[0.0,\"#0d0887\"],[0.1111111111111111,\"#46039f\"],[0.2222222222222222,\"#7201a8\"],[0.3333333333333333,\"#9c179e\"],[0.4444444444444444,\"#bd3786\"],[0.5555555555555556,\"#d8576b\"],[0.6666666666666666,\"#ed7953\"],[0.7777777777777778,\"#fb9f3a\"],[0.8888888888888888,\"#fdca26\"],[1.0,\"#f0f921\"]]},\"colorway\":[\"#636efa\",\"#EF553B\",\"#00cc96\",\"#ab63fa\",\"#FFA15A\",\"#19d3f3\",\"#FF6692\",\"#B6E880\",\"#FF97FF\",\"#FECB52\"],\"font\":{\"color\":\"#2a3f5f\"},\"geo\":{\"bgcolor\":\"white\",\"lakecolor\":\"white\",\"landcolor\":\"#E5ECF6\",\"showlakes\":true,\"showland\":true,\"subunitcolor\":\"white\"},\"hoverlabel\":{\"align\":\"left\"},\"hovermode\":\"closest\",\"mapbox\":{\"style\":\"light\"},\"paper_bgcolor\":\"white\",\"plot_bgcolor\":\"#E5ECF6\",\"polar\":{\"angularaxis\":{\"gridcolor\":\"white\",\"linecolor\":\"white\",\"ticks\":\"\"},\"bgcolor\":\"#E5ECF6\",\"radialaxis\":{\"gridcolor\":\"white\",\"linecolor\":\"white\",\"ticks\":\"\"}},\"scene\":{\"xaxis\":{\"backgroundcolor\":\"#E5ECF6\",\"gridcolor\":\"white\",\"gridwidth\":2,\"linecolor\":\"white\",\"showbackground\":true,\"ticks\":\"\",\"zerolinecolor\":\"white\"},\"yaxis\":{\"backgroundcolor\":\"#E5ECF6\",\"gridcolor\":\"white\",\"gridwidth\":2,\"linecolor\":\"white\",\"showbackground\":true,\"ticks\":\"\",\"zerolinecolor\":\"white\"},\"zaxis\":{\"backgroundcolor\":\"#E5ECF6\",\"gridcolor\":\"white\",\"gridwidth\":2,\"linecolor\":\"white\",\"showbackground\":true,\"ticks\":\"\",\"zerolinecolor\":\"white\"}},\"shapedefaults\":{\"line\":{\"color\":\"#2a3f5f\"}},\"ternary\":{\"aaxis\":{\"gridcolor\":\"white\",\"linecolor\":\"white\",\"ticks\":\"\"},\"baxis\":{\"gridcolor\":\"white\",\"linecolor\":\"white\",\"ticks\":\"\"},\"bgcolor\":\"#E5ECF6\",\"caxis\":{\"gridcolor\":\"white\",\"linecolor\":\"white\",\"ticks\":\"\"}},\"title\":{\"x\":0.05},\"xaxis\":{\"automargin\":true,\"gridcolor\":\"white\",\"linecolor\":\"white\",\"ticks\":\"\",\"title\":{\"standoff\":15},\"zerolinecolor\":\"white\",\"zerolinewidth\":2},\"yaxis\":{\"automargin\":true,\"gridcolor\":\"white\",\"linecolor\":\"white\",\"ticks\":\"\",\"title\":{\"standoff\":15},\"zerolinecolor\":\"white\",\"zerolinewidth\":2}}}},                        {\"responsive\": true}                    ).then(function(){\n",
       "                            \n",
       "var gd = document.getElementById('5b0d60c3-9682-4d91-a2dd-c0e82e7fcd49');\n",
       "var x = new MutationObserver(function (mutations, observer) {{\n",
       "        var display = window.getComputedStyle(gd).display;\n",
       "        if (!display || display === 'none') {{\n",
       "            console.log([gd, 'removed!']);\n",
       "            Plotly.purge(gd);\n",
       "            observer.disconnect();\n",
       "        }}\n",
       "}});\n",
       "\n",
       "// Listen for the removal of the full notebook cells\n",
       "var notebookContainer = gd.closest('#notebook-container');\n",
       "if (notebookContainer) {{\n",
       "    x.observe(notebookContainer, {childList: true});\n",
       "}}\n",
       "\n",
       "// Listen for the clearing of the current output cell\n",
       "var outputEl = gd.closest('.output');\n",
       "if (outputEl) {{\n",
       "    x.observe(outputEl, {childList: true});\n",
       "}}\n",
       "\n",
       "                        })                };                });            </script>        </div>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "### import plotly.graph_objects as go\n",
    "import numpy as np\n",
    "\n",
    "#size = 1/((quality_mesh.flatten()/np.max(quality_mesh.flatten()))**1)\n",
    "size = 1/quality_mesh.flatten()**5*5\n",
    "\n",
    "fig = go.Figure(data=[go.Scatter3d(x=rho_c_mesh.flatten(), \n",
    "                                   y=nu_mesh.flatten(), \n",
    "                                   z=zeta_mesh.flatten(),\n",
    "                                   mode='markers',\n",
    "                                  marker=dict(\n",
    "                                            size=size,\n",
    "                                            color=size,          # set color to an array/list of desired values\n",
    "                                    colorscale='inferno',   # choose a colorscale\n",
    "                                    opacity=0.8,\n",
    "                                      \n",
    "                                ))])\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1a7f17b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "0b2f2cf9",
   "metadata": {},
   "source": [
    "i had used the scaling func from Mendez-santos: where $x = np.power(l, -a) * y$ and $y = np.power(l, 1./nu) * (rho-rho_c)$. Now having changed to $x = np.power(l, -zeta/nu) * y$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 405,
   "id": "c2968441",
   "metadata": {},
   "outputs": [],
   "source": [
    "def goal_function_min(x):\n",
    "    rho_c = x[0] \n",
    "    nu = x[1]\n",
    "    zeta=x[2]\n",
    "\n",
    "    #XX = x_dat_2_new(rho=rho_mesh_local, rho_c=rho_c, l=l_mesh_local ,nu=nu) \n",
    "    XX = x_BKT(rho=rho_mesh_local, rho_c=rho_c, l=l_mesh_local ,nu=nu) \n",
    "    \n",
    "    \n",
    "    YY = y_dat_new(y=Y, l=l_mesh_local, zeta=zeta, nu=nu)\n",
    "    #print('y',YY)\n",
    "    dYY = y_dat_new(y=dY, l=l_mesh_local, zeta=zeta, nu=nu)\n",
    "    #print('dy',dYY)\n",
    "    \n",
    "    # sorting for BKT\n",
    "    sorting_index = np.argsort(XX)[0]\n",
    "    #print(sorting_index)\n",
    "    \n",
    "    XX = XX[:,sorting_index]\n",
    "    YY = YY[:,sorting_index]\n",
    "    dYY = dYY[:, sorting_index]\n",
    "    \n",
    "    \n",
    "    return quality(XX,YY,dYY)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 490,
   "id": "c5f5fd64",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "      fun: 0.9095750034698792\n",
       " hess_inv: <3x3 LbfgsInvHessProduct with dtype=float64>\n",
       "      jac: array([ 8.73819141, -0.11017891,  0.11013496])\n",
       "  message: 'ABNORMAL_TERMINATION_IN_LNSRCH'\n",
       "     nfev: 208\n",
       "      nit: 1\n",
       "     njev: 52\n",
       "   status: 2\n",
       "  success: False\n",
       "        x: array([3.5000079 , 2.00000198, 1.99999802])"
      ]
     },
     "execution_count": 490,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "minimize(goal_function_min,x0=[3.5,2.0,2.0], bounds=((3.3,3.9),(1.9,2.1),(1.9,2.1)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 470,
   "id": "8fb7d010",
   "metadata": {},
   "outputs": [],
   "source": [
    "import fssa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 471,
   "id": "ab606388",
   "metadata": {},
   "outputs": [],
   "source": [
    "def scaledata(l, rho, a, da, rho_c, nu, zeta):\n",
    "    r'''\n",
    "    Scale experimental data according to critical exponents\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    l, rho : 1-D array_like\n",
    "       finite system sizes `l` and parameter values `rho`\n",
    "\n",
    "    a, da : 2-D array_like of shape (`l`.size, `rho`.size)\n",
    "       experimental data `a` with standard errors `da` obtained at finite\n",
    "       system sizes `l` and parameter values `rho`, with\n",
    "       ``a.shape == da.shape == (l.size, rho.size)``\n",
    "\n",
    "    rho_c : float in range [rho.min(), rho.max()]\n",
    "       (assumed) critical parameter value with ``rho_c >= rho.min() and rho_c\n",
    "       <= rho.max()``\n",
    "\n",
    "    nu, zeta : float\n",
    "       (assumed) critical exponents\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    :py:class:`ScaledData`\n",
    "       scaled data `x`, `y` with standard errors `dy`\n",
    "\n",
    "    x, y, dy : ndarray\n",
    "       two-dimensional arrays of shape ``(l.size, rho.size)``\n",
    "\n",
    "    Notes\n",
    "    -----\n",
    "    Scale data points :math:`(\\varrho_j, a_{ij}, da_{ij})` observed at finite\n",
    "    system sizes :math:`L_i` and parameter values :math:`\\varrho_i` according\n",
    "    to the finite-size scaling ansatz\n",
    "\n",
    "    .. math::\n",
    "\n",
    "       L^{-\\zeta/\\nu} a_{ij} = \\tilde{f}\\left( L^{1/\\nu} (\\varrho_j -\n",
    "       \\varrho_c) \\right).\n",
    "\n",
    "    The output is the scaled data points :math:`(x_{ij}, y_{ij}, dy_{ij})` with\n",
    "\n",
    "    .. math::\n",
    "\n",
    "       x_{ij} & = L_i^{1/\\nu} (\\varrho_j - \\varrho_c) \\\\\n",
    "       y_{ij} & = L_i^{-\\zeta/\\nu} a_{ij} \\\\\n",
    "       dy_{ij} & = L_i^{-\\zeta/\\nu} da_{ij}\n",
    "\n",
    "    such that all data points :ref:`collapse <data-collapse-method>` onto the\n",
    "    single curve :math:`\\tilde{f}(x)` with the right choice of\n",
    "    :math:`\\varrho_c, \\nu, \\zeta` [4]_ [5]_.\n",
    "\n",
    "    Raises\n",
    "    ------\n",
    "    ValueError\n",
    "       If `l` or `rho` is not 1-D array_like, if `a` or `da` is not 2-D\n",
    "       array_like, if the shape of `a` or `da` differs from ``(l.size,\n",
    "       rho.size)``\n",
    "\n",
    "    References\n",
    "    ----------\n",
    "\n",
    "    .. [4] M. E. J. Newman and G. T. Barkema, Monte Carlo Methods in\n",
    "       Statistical Physics (Oxford University Press, 1999)\n",
    "\n",
    "    .. [5] K. Binder and D. W. Heermann, `Monte Carlo Simulation in Statistical\n",
    "       Physics <http://dx.doi.org/10.1007/978-3-642-03163-2>`_ (Springer,\n",
    "       Berlin, Heidelberg, 2010)\n",
    "    '''\n",
    "\n",
    "    # l should be 1-D array_like\n",
    "    l = np.asanyarray(l)\n",
    "    if l.ndim != 1:\n",
    "        raise ValueError(\"l should be 1-D array_like\")\n",
    "\n",
    "    # rho should be 1-D array_like\n",
    "    rho = np.asanyarray(rho)\n",
    "    if rho.ndim != 1:\n",
    "        raise ValueError(\"rho should be 1-D array_like\")\n",
    "\n",
    "    # a should be 2-D array_like\n",
    "    a = np.asanyarray(a)\n",
    "    if a.ndim != 2:\n",
    "        raise ValueError(\"a should be 2-D array_like\")\n",
    "\n",
    "    # a should have shape (l.size, rho.size)\n",
    "    if a.shape != (l.size, rho.size):\n",
    "        raise ValueError(\"a should have shape (l.size, rho.size)\")\n",
    "\n",
    "    # da should be 2-D array_like\n",
    "    da = np.asanyarray(da)\n",
    "    if da.ndim != 2:\n",
    "        raise ValueError(\"da should be 2-D array_like\")\n",
    "\n",
    "    # da should have shape (l.size, rho.size)\n",
    "    if da.shape != (l.size, rho.size):\n",
    "        raise ValueError(\"da should have shape (l.size, rho.size)\")\n",
    "\n",
    "    # rho_c should be float\n",
    "    rho_c = float(rho_c)\n",
    "\n",
    "    # rho_c should be in range\n",
    "    if rho_c > rho.max() or rho_c < rho.min():\n",
    "        warnings.warn(\"rho_c is out of range\", RuntimeWarning)\n",
    "\n",
    "    # nu should be float\n",
    "    nu = float(nu)\n",
    "\n",
    "    # zeta should be float\n",
    "    zeta = float(zeta)\n",
    "\n",
    "    l_mesh, rho_mesh = np.meshgrid(l, rho, indexing='ij')\n",
    "\n",
    "    x = np.power(l_mesh, 1. / nu) * (rho_mesh - rho_c)\n",
    "    x = l_mesh*np.exp(-1.4/np.sqrt(abs(rho_mesh-rho_c)))\n",
    "    y = np.power(l_mesh, - zeta / nu) * a\n",
    "    dy = np.power(l_mesh, - zeta / nu) * da\n",
    "\n",
    "    sorting_index = np.argsort(x)[0]\n",
    "    \n",
    "    x = x[:,sorting_index]\n",
    "    y = y[:,sorting_index]\n",
    "    dy = dy[:, sorting_index]\n",
    "\n",
    "    return ScaledData(x, y, dy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 472,
   "id": "84db8dfd",
   "metadata": {},
   "outputs": [
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-472-051f4fb10445>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mrho0\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinspace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m20\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m     \u001b[0mres\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfssa\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mautoscale\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mLs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrho\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mY\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdY\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mrho0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      6\u001b[0m     \u001b[0max\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mscatter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrho0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mres\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'rho'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mc\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'b'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m     \u001b[0max\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mscatter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrho0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mres\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'drho'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mc\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'b'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/anaconda3/lib/python3.8/site-packages/fssa/fssa.py\u001b[0m in \u001b[0;36mautoscale\u001b[0;34m(l, rho, a, da, rho_c0, nu0, zeta0, x_bounds, **kwargs)\u001b[0m\n\u001b[1;32m    609\u001b[0m         )\n\u001b[1;32m    610\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 611\u001b[0;31m     ret = scipy.optimize.minimize(\n\u001b[0m\u001b[1;32m    612\u001b[0m         \u001b[0mgoal_function\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    613\u001b[0m         \u001b[0;34m[\u001b[0m\u001b[0mrho_c0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnu0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mzeta0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/anaconda3/lib/python3.8/site-packages/scipy/optimize/_minimize.py\u001b[0m in \u001b[0;36mminimize\u001b[0;34m(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)\u001b[0m\n\u001b[1;32m    595\u001b[0m         \u001b[0;31m# custom method should be able to accept whatever bounds/constraints\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    596\u001b[0m         \u001b[0;31m# are provided to it.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 597\u001b[0;31m         return method(fun, x0, args=args, jac=jac, hess=hess, hessp=hessp,\n\u001b[0m\u001b[1;32m    598\u001b[0m                       \u001b[0mbounds\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mbounds\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mconstraints\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mconstraints\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    599\u001b[0m                       callback=callback, **options)\n",
      "\u001b[0;32m/opt/anaconda3/lib/python3.8/site-packages/fssa/optimize.py\u001b[0m in \u001b[0;36m_minimize_neldermead\u001b[0;34m(func, x0, args, callback, xtol, ftol, maxiter, maxfev, disp, return_all, **unknown_options)\u001b[0m\n\u001b[1;32m     88\u001b[0m         \u001b[0mxbar\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnumpy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreduce\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msim\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0mN\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     89\u001b[0m         \u001b[0mxr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mrho\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mxbar\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mrho\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0msim\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 90\u001b[0;31m         \u001b[0mfxr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mxr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     91\u001b[0m         \u001b[0mdoshrink\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     92\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/anaconda3/lib/python3.8/site-packages/scipy/optimize/optimize.py\u001b[0m in \u001b[0;36mfunction_wrapper\u001b[0;34m(*wrapper_args)\u001b[0m\n\u001b[1;32m    459\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mfunction_wrapper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mwrapper_args\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    460\u001b[0m         \u001b[0mncalls\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 461\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mfunction\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mwrapper_args\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    462\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    463\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mncalls\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfunction_wrapper\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/anaconda3/lib/python3.8/site-packages/fssa/fssa.py\u001b[0m in \u001b[0;36mgoal_function\u001b[0;34m(x)\u001b[0m\n\u001b[1;32m    605\u001b[0m             \u001b[0mrho\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mrho\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ml\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0ml\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mda\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mda\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnu\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mzeta\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrho_c\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    606\u001b[0m         )\n\u001b[0;32m--> 607\u001b[0;31m         return quality(\n\u001b[0m\u001b[1;32m    608\u001b[0m             \u001b[0mmy_x\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmy_y\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmy_dy\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx_bounds\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mx_bounds\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    609\u001b[0m         )\n",
      "\u001b[0;32m/opt/anaconda3/lib/python3.8/site-packages/fssa/fssa.py\u001b[0m in \u001b[0;36mquality\u001b[0;34m(x, y, dy, x_bounds)\u001b[0m\n\u001b[1;32m    433\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    434\u001b[0m             \u001b[0;31m# master curve estimate\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 435\u001b[0;31m             master_y[i, j], master_dy2[i, j] = _wls_linearfit_predict(\n\u001b[0m\u001b[1;32m    436\u001b[0m                 \u001b[0mx\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mj\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mw\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mw\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mwx\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mwx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mwy\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mwy\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mwxx\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mwxx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mwxy\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mwxy\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mselect\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mselect\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    437\u001b[0m             )\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsEAAAD8CAYAAACWynkuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAklUlEQVR4nO3df5Bdd3nf8ffHMoYIKCiV2nFtSzLEYFxIMGwNLUxDCBjhP2w6YVqZdTCMW00IpqTpdEqrmThjxjOhmRImEydhBzwmnS2GEJoqqTOuG+y6DcjxOhjjHzVRjCXLZWqBDUmq1MT20z/OEb5ar7R3d++Pc33er5k7957vOffej3b30X327Dnnm6pCkiRJ6pNTph1AkiRJmjSbYEmSJPWOTbAkSZJ6xyZYkiRJvWMTLEmSpN6xCZYkSVLv2ARLHZPkuiSPJrnnBOuT5FeTHEhyd5LXDay7PMmftrfLJ5da6ifrVZpdNsFS91wP7DrJ+ncC57S3PcBvACT5YeAq4A3ABcBVSbaMNamk67FepZlkEyx1TFXdBjx2kk0uAX6rGvuBlyY5HXgHcHNVPVZVjwM3c/IPZ0kbZL1Ks+vUcbzo1q1ba+fOneN4aWkm3Xnnnd+uqm0jerkzgIcHlg+3Yycaf5Yke2j2SvHCF77w9eeee+6Iokmzz3qVZsdG6nUsTfDOnTtZWloax0tLMynJwWlnGFRVC8ACwNzcXFmv0jOsV2l2bKRePRxCM2NxEXbuhFNOae4XF6edaGoeAc4aWD6zHTvRuKTpsV6ljrIJ1kxYXIQ9e+DgQahq7vfs6W0jvA94b3vW+RuB71XVt4CbgAuTbGlPsLmwHZM0Pdar1FFjORxCGrW9e+Ho0ePHjh5txufnp5NpXJJ8FngLsDXJYZozyJ8HUFW/CdwIXAQcAI4C72/XPZbko8Ad7UtdXVUnO2FH0gZZr9LsWrUJTvIC4Dbg+e32X6iqq8YdTBp06NDaxmdZVV26yvoCPniCddcB140jl6Rns16l2TXM4RBPAG+tqh8DXgvsav+kI03M9u1rG5ckSTqZVZvg9tqGf9kuPq+91VhTSctccw1s3nz82ObNzbgkSdJaDXViXJJNSe4CHqW5uPftK2yzJ8lSkqUjR46MOKb6bn4eFhZgxw5ImvuFhefe8cCSJGkyhmqCq+qpqnotzSVcLkjy6hW2Waiquaqa27ZtVNcYl54xPw8PPQRPP93c2wBLkqT1WtMl0qrqu8AtOLWjJEmSZtiqTXCSbUle2j7+IeDtwP8acy5JkiRpbIa5TvDpwGeSbKJpmj9fVb8/3liSJEnS+KzaBFfV3cD5E8giSZIkTYTTJkuSJKl3bIIlSZLUOzbBkiRJ6h2bYEmSJPWOTbC0TouLsHMnnHJKc7+4OO1EkiRpWMNcIk3SMouLsGcPHD3aLB882CyDM9lJkjQL3BMsrcPevc80wMccPdqMS5Kk7rMJltbh0KG1jUuSpG6xCZbWYfv2tY2vRZJdSR5IciDJR1ZY/ytJ7mpv30jy3YF1Tw2s27fxNJJWY81Ks8ljgqV1uOaa448JBti8uRnfiHZ68muBtwOHgTuS7Kuq+45tU1X/YmD7D3H8jI5/VVWv3VgKScOyZqXZ5Z5gaR3m52FhAXbsgKS5X1gYyUlxFwAHqurBqvo+cANwyUm2vxT47IbfVdJ6WbPSjLIJltZpfh4eegiefrq5H9FVIc4AHh5YPtyOPUuSHcDZwJcGhl+QZCnJ/iTvOtGbJNnTbrd05MiREcSWemvsNWu9SuNhEyzNrt3AF6rqqYGxHVU1B7wH+ESSl6/0xKpaqKq5qprbtm3bJLJKWmfNWq/SeNgES93yCHDWwPKZ7dhKdrPsz6pV9Uh7/yBwK8cfeyhp9KxZaUbZBEvdcgdwTpKzk5xG86H5rDPGk5wLbAG+MjC2Jcnz28dbgTcB9y1/rqSRsmalGeXVIaQOqaonk1wJ3ARsAq6rqnuTXA0sVdWxD9fdwA1VVQNPfxXwySRP0/yC+0uDZ6hLGj1rVppdNsFSx1TVjcCNy8Z+YdnyL67wvC8DrxlrOEnPYs1Ks8nDISRJktQ7NsGSJEnqHZtgSZIk9Y5NsCRJknrHJliSJEm9YxMsSZKk3rEJliRJUu/YBEuSJKl3bIIlSZLUOzbBkiRJ6h2bYEmSJPWOTbAkSZJ6xyZYkiRJvWMTLEmSpN6xCZYkSVLv2ARLkiSpd2yCpY5JsivJA0kOJPnICuvfl+RIkrva2z8dWHd5kj9tb5dPNrnUT9asNJtOnXYASc9Isgm4Fng7cBi4I8m+qrpv2aafq6orlz33h4GrgDmggDvb5z4+gehSL1mz0uxyT7DULRcAB6rqwar6PnADcMmQz30HcHNVPdZ+iN4M7BpTTkkNa1aaUas2wUnOSnJLkvuS3Jvkw5MIJvXUGcDDA8uH27HlfirJ3Um+kOSsNT6XJHuSLCVZOnLkyChyS3019pq1XqXxGGZP8JPAv6yq84A3Ah9Mct54Y0k6id8DdlbVj9LsOfrMWl+gqhaqaq6q5rZt2zbygJKOs6GatV6l8Vi1Ca6qb1XVn7SP/wK4nxPsXZK0YY8AZw0sn9mO/UBVfaeqnmgXPwW8ftjnSho5a1aaUWs6JjjJTuB84PYV1vnnGmnj7gDOSXJ2ktOA3cC+wQ2SnD6weDHNL6YANwEXJtmSZAtwYTsmaXysWWlGDX11iCQvAn4H+Lmq+vPl66tqAVgAmJubq5EllHqkqp5MciXNB+Em4LqqujfJ1cBSVe0D/nmSi2kOVXoMeF/73MeSfJTmQxng6qp6bOL/CKlHrFlpdqVq9X41yfOA3wduqqqPr7b93NxcLS0tjSCe9NyQ5M6qmpt2jpVYr9LxrFdpdmykXoe5OkSATwP3D9MAS5IkSV03zDHBbwJ+GnjrwGw3F405lyRJkjQ2qx4TXFX/E8gEskiSJEkT4YxxkiRJ6h2bYEmSJPWOTbAkSZJ6xyZYkiRJvWMTLEmSpN6xCdYJLS7Czp1wyinN/eLitBNJkiSNxtDTJqtfFhdhzx44erRZPniwWQaYn59eLkmSpFFwT7BWtHfvMw3wMUePNuOSJEmzziZYKzp0aG3jkiRJs8QmWCvavn1t45IkSbPEJlgruuYa2Lz5+LHNm5txSZKkWWcTrBXNz8PCAuzYAUlzv7DgSXGSJOm5wSZYJzQ/Dw89BE8/3dzbAE9Gkl1JHkhyIMlHVlj/80nuS3J3kj9MsmNg3VNJ7mpv+yabXOona1aaTV4iTeqQJJuAa4G3A4eBO5Lsq6r7Bjb7KjBXVUeTfAD4d8A/adf9VVW9dpKZpT6zZqXZ5Z7gDnFyCgEXAAeq6sGq+j5wA3DJ4AZVdUtVHbuA3X7gzAlnlPQMa1aaUTbBHXFscoqDB6HqmckpbIR75wzg4YHlw+3YiVwB/MHA8guSLCXZn+RdJ3pSkj3tdktHjhzZUGCp58Zes9arNB42wR3h5BRaqySXAXPALw8M76iqOeA9wCeSvHyl51bVQlXNVdXctm3bJpBW0npr1nqVxsMmuCOcnEKtR4CzBpbPbMeOk+RtwF7g4qp64th4VT3S3j8I3AqcP86wkqxZaVbZBHeEk1OodQdwTpKzk5wG7AaOO2M8yfnAJ2k+TB8dGN+S5Pnt463Am4DBk3MkjZ41K80om+COcHIKAVTVk8CVwE3A/cDnq+reJFcnubjd7JeBFwG/veyySq8ClpJ8DbgF+KVlZ6hLGjFrVppdXiKtI45dg3fv3uYQiO3bmwbYa/P2T1XdCNy4bOwXBh6/7QTP+zLwmvGmk7ScNSvNJpvgDpmft+mVJEmaBA+HkCRJUu/0ugl2cgpJkqR+6u3hEMcmpzh2bd5jk1OAhyRIkiQ91/V2T7CTU0iSJPVXb5tgJ6eQJEnqr942wU5OIUmS1F+9bYKdnEKSJKm/etsEz8/DwgLs2AFJc7+w4ElxkiRJfdDbq0OAk1NIkiT1VW/3BEuSJKm/bIIlSZLUOxNvgp2lTZIkSdM20WOCnaVNkiRJXTDRPcHO0iZJkqQumGgT7CxtkiRJ6oJVm+Ak1yV5NMk9G30zZ2mThpNkV5IHkhxI8pEV1j8/yefa9bcn2Tmw7t+04w8kecdEg0s9ZL1Ks2mYPcHXA7tG8WbO0iatLskm4FrgncB5wKVJzlu22RXA41X1I8CvAB9rn3sesBv4uzR1++vt60kaA+tVml2rNsFVdRvw2CjezFnapKFcAByoqger6vvADcAly7a5BPhM+/gLwE8mSTt+Q1U9UVXfBA60rydpPKxXaUaN7OoQSfYA7bUe+MskD6zylK0HD/Ltyy6Dyy4bVYp12Qp8e6oJGuY43nMtx441bHsG8PDA8mHgDSfapqqeTPI94G+24/uXPfeM5W+wrF6fGMXhTmPSlZ+DE+lyPrOt3yvXsK31erwuf2/Ntn5dzreWej3OyJrgqloAFobdPslSVc2N6v3XyxzmmIUcozZYr13+N3Y5G3Q7n9nWL8nStDMMmpV6hW7nM9v6dTnfRurVGeOk7nkEOGtg+cx2bMVtkpwKvAT4zpDPlTQ61qs0o2yCpe65AzgnydlJTqM5cWbfsm32AZe3j98NfKmqqh3f3Z6NfjZwDvDHE8ot9ZH1Ks2oVQ+HSPJZ4C3A1iSHgauq6tMjeO+hD50YM3MczxzHm3iO9pjBK4GbgE3AdVV1b5KrgaWq2gd8GvgPSQ7QnLi6u33uvUk+D9wHPAl8sKqeWuUtu/K1XkmXs0G385lt/YbOZ70+S5fzmW39upxv3dnS/DIqSZIk9YeHQ0iSJKl3bIIlSZLUO2NtglebcjmNX22njLw7yeumlGO+ff+vJ/lykh+bRo6B7f5ekieTvHtaOZK8JcldSe5N8t+nkSPJS5L8XpKvtTneP4YMZyW5Jcl97Xt8eIVtJvJzOm4bmdq1A9l+vv0e3Z3kD5Os5brLY802sN1PJakkE72M0DD5kvzjgZ/x/9iVbEm2t/X31fZ7e9EEs3Xi8+kk+azXMeUb2G7iNdvleh0m37Rqdmz1WlVjuwH/EHgdcM8J1l8E/AEQ4I3A7VPK8Q+ALe3jd04rR7vNJuBLwI3Au6f09XgpzYka29vlvzWlHP8W+Fj7eBvNCSWnjTjD6cDr2scvBr4BnLdsm4n8nI7z1v5c/RnwMuA04Gsr/Dt/FvjN9vFu4HMdyvYTwOb28Qe6lG3gZ+c2mokP5jr2fT0H+OrA/3Fjqed1ZlsAPtA+Pg94aIJfu058Pm3ga2e9rjNfu93Ea7bL9bqGfFOp2XHV61j3BNfqUy5fAvxWNfYDL01y+qRzVNWXq+rxdnE/zbUaR26IrwfAh4DfAR4dR4Yhc7wH+GJVHWq3H0uWIXIU8OIkAV7UbvvkiDN8q6r+pH38F8D9PHvGpon8nI7ZRqZ2nXq2qrqlqo62i2Or0fVka30U+Bjw/yaU65hh8v0z4Npj/8eNq57Xma2Av9E+fgnwvyeUrTOfTydgvY4xX2saNdvleh0231Rqdlz1Ou1jgleabvJZU0ZO2BU0v01MXJIzgH8E/MY03n/AK4AtSW5NcmeS904px68Br6Ipsq8DH66qp8f1Zu2fE88Hbl+2qos/p2s1zL/huKldgWNTu3Yh26BJ1uiq2do/u51VVf9lQpkGDfO1ewXwiiR/lGR/kl0dyvaLwGVpLr95I81OgK6YZt1br+vX5Zrtcr3CbNfsuup1ZNMmPxck+Qmagn3zlCJ8AvjXVfX0ZH6hP6FTgdcDPwn8EPCVJPur6hsTzvEO4C7grcDLgZuT/I+q+vNRv1GSF9Hsgf+5cby+RiPJZcAc8OPTzgKQ5BTg48D7phzlZE6l+RPrW2j2yN2W5DVV9d1phmpdClxfVf8+yd+nuZbuq8f5y64mp2v1CjNRs12uV3iO1ey09wR3ZsrIJD8KfAq4pKq+M40MNP9Z3JDkIZpZhX49ybumkOMwcFNV/d+q+jbNcVNjOVlwFe+nOSyjquoA8E3g3FG/SZLn0TTAi1X1xRU26czP6QZsZGrXLmQjyduAvcDFVfXEBHINk+3FwKuBW9u6fSOwb4In2gzztTsM7Kuqv66qb9Ic935OR7JdAXweoKq+ArwA2DqBbMOYZt1br+vX5Zrtcr0Om6+rNbuuep12E7wPeG97Vt8bge9V1bcmHSLJduCLwE9PYW/nD1TV2VW1s6p20hzj9bNV9btTiPKfgTcnOTXJZuANNMfKTtohmr3RJPnbwCuBB0f5Bu0xdJ8G7q+qj59gs078nG7QRqZ2nXq2JOcDn6T5QJ3kMXInzVZV36uqrQN1u7/NuNSFfK3fpdmrRJKtNH9uHWkdbSDbYI2/iuYD9cgEsg1jmnVvvY4p35Rrtsv1Omy+rtbs+up1+Zlyo7wBnwW+Bfw1zW83VwA/A/xMuz7AtTRnI36dMZ2hOUSOTwGP0/zp/S6aqS4nnmPZttczvqtDrJoD+Fc0V4i4h+YQgWl8X/4O8F/bn417gMvGkOHNNAf63z3w/b9oGj+n4761/65vtP+Ove3Y1TQfAND8Z/bbwAHgj4GXdSjbfwP+z8D3aF9Xsi3b9tZJ/3wM8bULzZ9/72t/fnd3KNt5wB/RnIV+F3DhBLN14vNpA18763Wd+ZZtO9Ga7XK9DplvKjU7rnp12mRJkiT1zrQPh5C0zEYuCp7k8iR/2t4uX+n5kkbHepVml02w1D3XAye7LM47aU6UOAfYQ3tJvSQ/DFxFcwz3BcBVSbaMNamk67FepZlkEyx1TK3/ouDvAG6uqsequdD6zZz8w1nSBlmv0uway3WCt27dWjt37hzHS0sz6c477/x2VW0b0cud6KLgQ18sPMkemr1SvPCFL3z9ueeO/Mpz0syyXqXZsZF6XbUJTvICmuvEPr/d/gtVddXJnrNz506WliZ1hSAtt7gIe/fCoUOwfTtccw3Mz087Vb8lOTjtDIOqaoFmDnjm5ubKepWeYb1Ks2Mj9TrM4RBPAG+tqh8DXgvsaq/Bpg5aXIQ9e+DgQahq7vfsacb1nHGii4I/Fyb1kJ5rrFepo1ZtgtvjmP6yXXxee/O6ah21dy8cPXr82NGjzbieM050UfCbgAuTbGlPsLmwHZM0Pdar1FFDHROcZBNwJ/AjwLVVdfsK2/zgmKXt27ePMqPW4NChtY2re5J8lmbGoK1JDtOcQf48gKr6TeBGmguaHwCO0kwvTVU9luSjNLP+AFxdVSc7YUfSBlmv0uwaqgmuqqeA1yZ5KfCfkry6qu5Zts1xxyyNOqiGs317cwjESuOaDVV16SrrC/jgCdZdB1w3jlySns16lWbXmi6RVlXfBW7By7h01jXXwObNx49t3tyMS5IkqbFqE5xkW7sHmCQ/BLwd+F9jzqV1mp+HhQXYsQOS5n5hwatDSJIkDRrmcIjTgc+0xwWfAny+qn5/vLG0EfPzNr2SJEkns2oTXFV3A+dPIIskSZI0EU6bLEmSpN6xCZYkSVLv2ARLkiSpd2yCJUmS1Ds2wZIkSeodm2BJkiT1jk2wJEmSescmWJIkSb1jEyxJkqTesQmWJElS79gES5IkqXdsgiVJktQ7NsFSxyTZleSBJAeSfGSF9b+S5K729o0k3x1Y99TAun0TDS71lDUrzaZTpx1A0jOSbAKuBd4OHAbuSLKvqu47tk1V/YuB7T8EnD/wEn9VVa+dUFyp96xZaXa5J1jqlguAA1X1YFV9H7gBuOQk218KfHYiySStxJqVZpRNsNQtZwAPDywfbseeJckO4GzgSwPDL0iylGR/kned6E2S7Gm3Wzpy5MgIYku9NfaatV6l8bAJltZpcRF27oRTTmnuFxcnHmE38IWqempgbEdVzQHvAT6R5OUrPbGqFqpqrqrmtm3bNomsktZZs9arNB42wdI6LC7Cnj1w8CBUNfd79oykEX4EOGtg+cx2bCW7WfZn1ap6pL1/ELiV4489lDR61qw0o2yCpXXYuxeOHj1+7OjRZnyD7gDOSXJ2ktNoPjSfdcZ4knOBLcBXBsa2JHl++3gr8CbgvuXPlTRS1qw0o7w6hLQOhw6tbXxYVfVkkiuBm4BNwHVVdW+Sq4Glqjr24bobuKGqauDprwI+meRpml9wf2nwDHVJo2fNSrPLJlhah+3bm0MgVhrfqKq6Ebhx2dgvLFv+xRWe92XgNRtPIGktrFlpNnk4hLQO11wDmzcfP7Z5czMuSZK6zyZYWof5eVhYgB07IGnuFxaacUmS1H0eDiGt0/y8Ta8kSbPKPcGSJEnqnVWb4CRnJbklyX1J7k3y4UkEkyRJksZlmMMhngT+ZVX9SZIXA3cmudnLuEiSJGlWrbonuKq+VVV/0j7+C+B+TjAvuiRJkjQL1nRMcJKdNFM63r7Cuj1JlpIsHTlyZETxJEmSpNEbuglO8iLgd4Cfq6o/X76+qhaqaq6q5rZt2zbKjJIkSdJIDdUEJ3keTQO8WFVfHG8kSZIkabyGuTpEgE8D91fVx8cfSZIkSRqvYfYEvwn4aeCtSe5qbxeNOZckSZI0NqteIq2q/ieQCWSRJEmSJsIZ4yRJktQ7NsGSJEnqHZtgqWOS7EryQJIDST6ywvr3JTkycIz+Px1Yd3mSP21vl082udRP1qw0m4aZNlnShCTZBFwLvB04DNyRZN8K05R/rqquXPbcHwauAuaAopnifF9VPT6B6FIvWbPS7HJPsNQtFwAHqurBqvo+cANwyZDPfQdwc1U91n6I3gzsGlNOSQ1rVppRNsFSt5wBPDywfLgdW+6nktyd5AtJzlrjc53mXBqdsdes9SqNh02wNHt+D9hZVT9Ks+foM2t9Aac5lyZqQzVrvUrjYRMsdcsjwFkDy2e2Yz9QVd+pqifaxU8Brx/2uZJGzpqVZpRNsNQtdwDnJDk7yWnAbmDf4AZJTh9YvBi4v318E3Bhki1JtgAXtmOSxsealWaUV4eQOqSqnkxyJc0H4Sbguqq6N8nVwFJV7QP+eZKLgSeBx4D3tc99LMlHaT6UAa6uqscm/o+QesSalWZXqmrkLzo3N1dLS0sjf11pViW5s6rmpp1jJdardDzrVZodG6lXD4eQJElS79gES5IkqXdsgiVJktQ7NsGSJEnqHZtgSZIk9Y5NsCRJknrHJliSJEm9YxMsSZKk3rEJliRJUu/YBEuSJKl3bIIlSZLUOzbBkiRJ6h2bYEmSJPWOTbAkSZJ6Z9UmOMl1SR5Ncs8kAkmSJEnjNsye4OuBXWPOIUmSJE3Mqk1wVd0GPDaBLJKAJLuSPJDkQJKPrLD+55Pcl+TuJH+YZMfAuqeS3NXe9k02udRP1qw0m04d1Qsl2QPsAdi+ffuoXlbqlSSbgGuBtwOHgTuS7Kuq+wY2+yowV1VHk3wA+HfAP2nX/VVVvXaSmaU+s2al2TWyE+OqaqGq5qpqbtu2baN6WalvLgAOVNWDVfV94AbgksENquqWqjraLu4HzpxwRknPsGalGeXVIaRuOQN4eGD5cDt2IlcAfzCw/IIkS0n2J3nXiZ6UZE+73dKRI0c2FFjqubHXrPUqjcfIDoeQNFlJLgPmgB8fGN5RVY8keRnwpSRfr6o/W/7cqloAFgDm5uZqIoGlnltvzVqv0ngMc4m0zwJfAV6Z5HCSK8YfS+qtR4CzBpbPbMeOk+RtwF7g4qp64th4VT3S3j8I3AqcP86wkqxZaVYNc3WIS6vq9Kp6XlWdWVWfnkQwqafuAM5JcnaS04DdwHFnjCc5H/gkzYfpowPjW5I8v328FXgTMHhyjqTRs2alGeXhEFKHVNWTSa4EbgI2AddV1b1JrgaWqmof8MvAi4DfTgJwqKouBl4FfDLJ0zS/4P7SsjPUJY2YNSvNLptgqWOq6kbgxmVjvzDw+G0neN6XgdeMN52k5axZaTZ5dQhJkiT1jk2wJEmSescmWJIkSb1jEyxJkqTesQmWJElS79gES5IkqXdsgiVJktQ7NsGSJEnqHZtgSZIk9Y5NsCRJknrHJliSJEm9M/EmeHERdu6EU05p7hcXJ51AkiRJfXfqJN9scRH27IGjR5vlgwebZYD5+UkmkSRJUp9NdE/w3r3PNMDHHD3ajEuSJEmTMtEm+NChtY1LkiRJ4zDRJnj79rWNS5IkSeMw0Sb4mmtg8+bjxzZvbsYlPSPJriQPJDmQ5CMrrH9+ks+1629PsnNg3b9pxx9I8o6JBpd6yHqVZtNEm+D5eVhYgB07IGnuFxY8KU4alGQTcC3wTuA84NIk5y3b7Arg8ar6EeBXgI+1zz0P2A38XWAX8Ovt60kaA+tVml0Tv0Ta/Dw89BA8/XRzbwMsPcsFwIGqerCqvg/cAFyybJtLgM+0j78A/GSStOM3VNUTVfVN4ED7epLGw3qVZtRYLpF25513fjvJwVU22wp8exzvP2KzkhNmJ+us5ITRZd2xhm3PAB4eWD4MvOFE21TVk0m+B/zNdnz/sueesfwNkuwB2gsU8kSSe9aQb5K6/rPS5XxmW79XrmFb6/V4Xf7emm39upxvLfV6nLE0wVW1bbVtkixV1dw43n+UZiUnzE7WWckJs5V1LapqAViAbv8bu5wNup3PbOuXZGnaGQbNSr1Ct/OZbf26nG8j9eq0yVL3PAKcNbB8Zju24jZJTgVeAnxnyOdKGh3rVZpRNsFS99wBnJPk7CSn0Zw4s2/ZNvuAy9vH7wa+VFXVju9uz0Y/GzgH+OMJ5Zb6yHqVZtREp01eZmGK770Ws5ITZifrrOSEKWRtjxm8ErgJ2ARcV1X3JrkaWKqqfcCngf+Q5ADwGM0HL+12nwfuA54EPlhVT63yll3+fnQ5G3Q7n9nWb+h81uuzdDmf2davy/nWnS3NL6OSJElSf3g4hCRJknrHJliSJEm9M9YmOMl1SR490TUN0/jVdsrIu5O8bpx5TmSInPNtvq8n+XKSH5t0xoEsJ806sN3fS/JkkndPKtuy9181Z5K3JLkryb1J/vsk8y3Lsdr3/yVJfi/J19qs7590xlHYyNSuHcj280nua+vwD5Os5brLY802sN1PJakkE72M0DD5kvzj9ut3b5L/2JVsSbYnuSXJV9vv7UUTzNbpzyfrdXz5BrabeM12uV6HyTetmh1bvVbV2G7APwReB9xzgvUXAX8ABHgjcPs482wg5z8AtrSP3zmtnMNkbbfZBHwJuBF4dxdzAi+lORlke7v8t7r6NQX+LfCx9vE2mhNbTptW3nX+GzcBfwa8DDgN+Bpw3rJtfhb4zfbxbuBzHcr2E8Dm9vEHupSt3e7FwG00Ex/Mdez7eg7w1YH/wyZSa0NmWwA+0D4+D3hogl+7zn4+Wa/jzdduN/Ga7XK9riHfVGp2XPU61j3BVXUbTcNwIpcAv1WN/cBLk5w+zkwrWS1nVX25qh5vF/fTXMtxKob4mgJ8CPgd4NHxJ1rZEDnfA3yxqg6123c5awEvThLgRe22T04i2whtZGrXqWerqluq6mi7OMkaHObrBvBR4GPA/5tQrmOGyffPgGuP/R82wVobJlsBf6N9/BLgf08oW9c/n6zXMeZrTaNmu1yvw+abSs2Oq16nfUzwStNNPmvKyI65gua3jU5Kcgbwj4DfmHaWVbwC2JLk1iR3JnnvtAOdxK8Br6Ip9q8DH66qp6cbac2GqbXjpnYFjk3t2oVsgyZZg6tma//sdlZV/ZcJZRo0zNfuFcArkvxRkv1JdnUo2y8ClyU5TPOXqw9NJtpQpvn5ZL2uX5drtsv1CrNds+uq12leJ3jmJPkJmoJ+87SznMQngH9dVU9PZqfAup0KvB74SeCHgK8k2V9V35hurBW9A7gLeCvwcuDmJP+jqv58qql6KMllwBzw49POApDkFODjwPumHOVkTqX5E+tbaPbI3ZbkNVX13WmGal0KXF9V/z7J36e5lu6rZ/CXTK2ga/UKM1GzXa5XeI7V7LT3BM/MlJFJfhT4FHBJVX1n2nlOYg64IclDNDMT/XqSd0010coOAzdV1f+tqm/THJs1tRMOV/F+mkM3qqoOAN8Ezp1yprXayNSuXchGkrcBe4GLq+qJCeQaJtuLgVcDt7Y190Zg3wRPtBnma3cY2FdVf11V3wS+QfMh24VsVwCfB6iqrwAvALZOINswpvn5ZL2uX5drtsv1Omy+rtbsuup12k3wPuC97Vl9bwS+V1XfmnKmZ0myHfgi8NMd3VP5A1V1dlXtrKqdNMeJ/WxV/e50U63oPwNvTnJqks3AG4D7p5zpRA7R7LEmyd8GXgk8ONVEa7eRqV2nni3J+cAnaT5QJ3mM3EmzVdX3qmrrQM3tbzMudSFf63dp9iqRZCvNn1sn8fM7TLbB2noVzQfqkQlkG8Y0P5+s1zHlm3LNdrleh83X1ZpdX70Oc/bcem/AZ4FvAX9N89vNFcDPAD/Trg9wLc3ZiF9ngmdVrzHnp4DHaf4kfhfNVJgTzzlM1mXbXs/0rg6xak7gX9FcIeIe4Oe6+jUF/g7wX9uf0XuAy6aVdYP/zoto9ir8GbC3Hbua5gMAmv/Mfhs4APwx8LIOZftvwP8ZqMF9Xcm2bNtbJ/3/2BBfu9D8+fe+9md4d4eynQf8Ec1Z6HcBF04wW6c/n6zX8eVbtu1Ea7bL9TpkvqnU7Ljq1WmTJUmS1DvTPhxCkiRJmjibYEmSJPWOTbAkSZJ6xyZYkiRJvWMTLEmSpN6xCZYkSVLv2ARLkiSpd/4/GYDsyzLbWU4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x288 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "fig, ax = plt.subplots(2,3, figsize=(12,4))\n",
    "\n",
    "\n",
    "for rho0 in np.linspace(1,5,20):\n",
    "    res = fssa.autoscale(Ls, rho, Y, dY,rho0,2,2)\n",
    "    ax[0,0].scatter(rho0, res['rho'], c='b',)\n",
    "    ax[1,0].scatter(rho0, res['drho'], c='b')\n",
    "    \n",
    "print('rho done')\n",
    "for nu0 in np.linspace(1,5,20):\n",
    "    res = fssa.autoscale(Ls, rho, Y, dY,3.6,nu0,2)\n",
    "    ax[0,1].scatter(nu0, res['nu'],c='r')\n",
    "    ax[1,1].scatter(nu0, res['dnu'],c='r')\n",
    "\n",
    "print('nu done')\n",
    "for zeta0 in np.linspace(1,5,20):\n",
    "    res = fssa.autoscale(Ls, rho, Y, dY,3.6,2,zeta0)\n",
    "    ax[0,2].scatter(zeta0, res['zeta'],c='g')\n",
    "    ax[1,2].scatter(zeta0, res['dzeta'],c='g')\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 491,
   "id": "c1ddfbdf",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/lib/python3.8/site-packages/fssa/fssa.py:505: RuntimeWarning:\n",
      "\n",
      "invalid value encountered in sqrt\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "           dnu: 1.2761443204687717\n",
       "          drho: nan\n",
       "         dzeta: 2.3209718565489075\n",
       "        errors: array([       nan, 1.27614432, 2.32097186])\n",
       " final_simplex: (array([[3.91993678, 1.7985039 , 1.7442203 ],\n",
       "       [3.9151196 , 1.80343237, 1.75199902],\n",
       "       [3.92285265, 1.80215994, 1.74215271],\n",
       "       [3.91943889, 1.80470511, 1.75072343]]), array([0.06940015, 0.06940788, 0.06941019, 0.06941343]))\n",
       "           fun: 0.06940015457037921\n",
       "       message: 'Optimization terminated successfully.'\n",
       "          nfev: 80\n",
       "           nit: 43\n",
       "            nu: 1.7985038970242853\n",
       "           rho: 3.9199367832344016\n",
       "        status: 0\n",
       "       success: True\n",
       "         varco: array([[-9.96435222, -4.76961027, -2.81822282],\n",
       "       [-4.76961027,  1.62854433,  3.47746285],\n",
       "       [-2.81822282,  3.47746285,  5.38691036]])\n",
       "             x: array([3.91993678, 1.7985039 , 1.7442203 ])\n",
       "          zeta: 1.74422030018006"
      ]
     },
     "execution_count": 491,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fssa.autoscale(Ls, rho, Y, dY,3.6,2,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0ffb45d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "c5159811",
   "metadata": {},
   "source": [
    "Anton should:\n",
    "\n",
    "    - [] More BKT scaling\n",
    "    - [] look if others use the same seeds for diff. disorders."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 500,
   "id": "e925a525",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1. , 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2. , 2.1, 2.2,\n",
       "       2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3. , 3.1, 3.2, 3.3, 3.4, 3.5,\n",
       "       3.6, 3.7, 3.8, 3.9, 4. , 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8,\n",
       "       4.9, 5. , 5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6. ])"
      ]
     },
     "execution_count": 500,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linspace(1,6.,51)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 504,
   "id": "5c03899f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([3.3 , 3.33, 3.36, 3.39, 3.42, 3.45, 3.48, 3.51, 3.54, 3.57, 3.6 ,\n",
       "       3.63, 3.66, 3.69, 3.72, 3.75, 3.78, 3.81, 3.84, 3.87, 3.9 ])"
      ]
     },
     "execution_count": 504,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linspace(3.3,3.9,21)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8cef7a1f",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}