{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "7ec1ce2d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import networkx as nx\n",
    "import numpy as np\n",
    "from collections import Counter\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "69f7088e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# make a random spin vector\n",
    "def random_state():\n",
    "    return 2*((np.random.random(N) > 0.5).astype('float') - 0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d0aa5b84",
   "metadata": {},
   "outputs": [],
   "source": [
    "# cost function of spin glass\n",
    "def cost(J, x):\n",
    "    return (x.T @ J @ x)/N**1.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "68de9506",
   "metadata": {},
   "outputs": [],
   "source": [
    "# find best local update\n",
    "def local_search(J, f, x, verbose=False):\n",
    "    pr = print if verbose else lambda *x: x\n",
    "    best = f(J, x)\n",
    "    loops = 0\n",
    "    pr('-')\n",
    "    pr('start cost:', best)\n",
    "    while True:\n",
    "        if loops % 20 == 0:\n",
    "            pr('loop', loops, 'cost:', best)\n",
    "        bestpos = None\n",
    "        for i in range(N):\n",
    "            x[i] *= -1\n",
    "            newcost = f(J, x)\n",
    "            if newcost < best:\n",
    "                best = newcost\n",
    "                bestpos = i\n",
    "            x[i] *= -1\n",
    "        if bestpos == None:\n",
    "            break\n",
    "        x[bestpos] *= -1\n",
    "        loops += 1\n",
    "    pr('final cost:', best)\n",
    "    return x, best"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e663a5b0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-\n",
      "start cost: 0.08759910087678596\n",
      "loop 0 cost: 0.08759910087678596\n",
      "loop 20 cost: -0.39087070518004546\n",
      "loop 40 cost: -0.6478017710308193\n",
      "loop 60 cost: -0.8036595333440774\n",
      "loop 80 cost: -0.8964781476205971\n",
      "loop 100 cost: -0.9635012906426891\n",
      "final cost: -0.9850169828722638\n"
     ]
    }
   ],
   "source": [
    "# initialize random couplings\n",
    "N = 200\n",
    "J = np.random.normal(size=(N,N))\n",
    "x = random_state()\n",
    "x, best = local_search(J, cost, x, verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b639fff",
   "metadata": {},
   "source": [
    "At large N, this should approach the Parisi constant $\\sqrt{2} * 0.7632... = 1.079...$, 2**0.5 * 0.7632but the exact value depends on the spin glass. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "259d1d4a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0793277908031462"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2**0.5 * 0.7632"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4683c7c7",
   "metadata": {},
   "source": [
    "Now we try using a bounded alternative function, $K' =  -tanh$. The strategy is to local search, then use tanh just below that area.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "de89e848",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7feee4933310>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAgj0lEQVR4nO3deXRd5X3u8e+j0aNky5aFLcsDYAKeAWGmDMwFGuLkJiFwW8pKm+tFCmmTld6GNrfDpXfgJm2mhhvqAAm9oaEkhEIKCYObMIQAlg0YgzE2NnjEludRtmX97h/nCA6yZEs+R9pneD5rnaW93/3uc34Y+Tze+917v4oIzMysdJUlXYCZmSXLQWBmVuIcBGZmJc5BYGZW4hwEZmYlriLpAo7H6NGjY9KkSUmXYWZWUBYtWrQlIuq7thdkEEyaNImWlpakyzAzKyiS3u6u3aeGzMxKnIPAzKzEOQjMzEqcg8DMrMQ5CMzMSlxOgkDSXZI2S1raw3ZJ+o6klZKWSDojY9vlkpant92ci3rMzKz3cnVE8EPg8qNsvwKYkn7NA74HIKkcuC29fSpwraSpOarJzMx6ISf3EUTEU5ImHaXLXOCfI/XM6+ckjZA0FpgErIyIVQCS7k33fS0XdXW1YNkmXl67470G6b3FI5sQ6q5rRut77crs0M323rxfb/q+/73VYz1Hth/jv6WbP4ueaioTVJaXUVEuqsrLqCgvo7LLcmV5GZXlZQyqLGNodQXDqiuorijr8c/JzJIzUDeUNQJrM9bXpdu6az+7uzeQNI/U0QQTJkw4riKefKOV//dc6n4KT8Mw8CrL9W4oDKuuoG5oFWOGVzOmZhD1w6oZU1PNhLohnFg/jNrBlUmXa1YyBioIuvtnYByl/cjGiPnAfIDm5ubj+hq/Ze50bpk7vcftmZP0ZAZF9NSn277dv8f7P+fI/n35vPf178VnH2/9PSxyuCNoPxwcPNxBe0cHh9rTy4c7OHQ4ONTRwaH21PL+Q4fZe6CdPenX3gPt7GlrZ1dbO9v2HmDRmu1s3nWAA+0dZBo9rIoTRw9jxvhaTp8wgtMnjGRc7SAfUZj1g4EKgnVAU8b6eGADUNVDeyJ6On3SpdeA1FJKIoJdbe1s3tXGW1v3sap1D6ta97Ji825+9Nzb3PnMagAmjhrCRaeO4ZLTGjh7ch0V5b7ozSwXBioIHgJuSo8BnA3sjIiNklqBKZImA+uBa4D/PEA1WZ6QRO3gSmoHVzKlYTjQ8O62Q4c7WLZxF4ve3s5Tb7Ryz/Nr+MFv3qKhpppPn9nEZ85qoqluSHLFmxUB5WLOYkk/Bi4ARgObgL8BKgEi4nal/qn9XVJXFu0DPhsRLel9rwS+BZQDd0XE/zzW5zU3N4cfOlea9h1s56k3WrmvZR2/Xr4ZSXx8diM3XXQyk0cPTbo8s7wmaVFENB/RXoiT1zsIDGDjzv3c8fRq7nn+bQ4dDv7w/El88ZJTGFpdkA/VNet3PQWBT7JawRpbO5i/+uhUnv7zi/j0meP5/tOrufQbT/LbN7cmXZpZQXEQWMGrH17NrZ+cyf2fP5dBleX83h3P8Z0FK+joKLyjXbMkOAisaJw5sY6ff+GDfGzWOL7x+Bt86b6XONjlslQzO5JPplpRGVpdwTc/M5spDcP5+qPL2bb3IN//g2YGVZYnXZpZ3vIRgRUdSdx44cl87VMzeWblFm68ZzGHDvvIwKwnDgIrWlc3N/F3c6ez4PXNfOWnSyjEK+TMBoJPDVlR+/1zJrJ1z0G++cQbTB1Xw+c+dGLSJZnlHR8RWNH7k4tP5vJpJ/C/f/E6z765JelyzPKOg8CKniT+/upZTBw1hC/f9zI79x9KuiSzvOIgsJIwrLqCb1w9m827D3DLz/tluguzguUgsJIxu2kEf3zBSdy/OPWcIjNLcRBYSfnCRVM4cfRQbvn5a77ZzCzNQWAlpaqijL+6aiqrtuzlh8+uTrocs7zgILCSc+EHxnDxqWP4zoKVbN97MOlyzBLnILCS9JUrTmXvwXa+//SqpEsxS5yDwErSKQ3DuWrmOH747Fts3XMg6XLMEpWTIJB0uaTlklZKurmb7f9V0kvp11JJhyXVpbe9JemV9DbPNmMD5k8unkLbocPMf8pHBVbasg4CSeXAbcAVwFTgWklTM/tExNcjYnZEzAb+AngyIrZldLkwvf2ImXPM+svJY4bxuzPH8S/Pr2HPgfakyzFLTC6OCOYAKyNiVUQcBO4F5h6l/7XAj3PwuWZZ+6MPTmb3gXbuW7g26VLMEpOLIGgEMv8WrUu3HUHSEFIT2N+f0RzAY5IWSZrX04dImiepRVJLa2trDso2S91kdubEkfzg2dUc9oxmVqJyEQTqpq2nv1FXAb/pclro/Ig4g9SppRslfbi7HSNifkQ0R0RzfX19dhWbZfijD05m7bb9PLFsU9KlmCUiF0GwDmjKWB8PbOih7zV0OS0UERvSPzcDD5A61WQ2YC6b2kBDTTX/6tNDVqJyEQQLgSmSJkuqIvVl/1DXTpJqgY8AD2a0DZU0vHMZuAxYmoOazHqtoryMT505nl8v38w7O9uSLsdswGUdBBHRDtwEPAosA+6LiFcl3SDphoyunwAei4i9GW0NwDOSXgZeAB6OiF9mW5NZX13d3ERHwE8X+ajASo8Kcfq+5ubmaGnxLQeWW9fM/y0bdrTx6z+7gLKy7oa+zAqbpEXdXabvO4vN0q5ubmLNtn0sXrM96VLMBpSDwCztsmknUF1Rxs9f7ulaB7Pi5CAwSxtWXcFFp47hkaXv+J4CKykOArMMH505jtbdB3h+9dakSzEbMA4CswwXnTqGIVXl/PzljUmXYjZgHARmGQZXlXPJaQ08+qpPD1npcBCYdXHZtAa27T3Ii756yEqEg8Csiw+fUk9FmXhi2eakSzEbEA4Csy5qBlVy9ol1fgidlQwHgVk3LjmtgZWb9/D21r3H7mxW4BwEZt245LQGAJ8espLgIDDrRlPdEE5pGMYCnx6yEuAgMOvBR06pp+Wt7ew/eDjpUsz6lYPArAcfnFLPwcMdvsvYip6DwKwHcybVUVVexjMrtiRdilm/chCY9WBwVTnNk0byzEoHgRW3nASBpMslLZe0UtLN3Wy/QNJOSS+lX3/d233NkvTBKaN5/Z3dbN7tKSyteGUdBJLKgduAK4CpwLWSpnbT9emImJ1+3dLHfc0S8aGT6wH4jY8KrIjl4ohgDrAyIlZFxEHgXmDuAOxr1u+mjath5JBKnvY4gRWxXARBI5A54/e6dFtX50p6WdIvJE3r475ImiepRVJLa2trDso2O7ayMnHuSaN4ftW2pEsx6ze5CILuZvnu+vzexcDEiJgF/CPwb33YN9UYMT8imiOiub6+/nhrNeuzOZPqWL9jP+u270u6FLN+kYsgWAc0ZayPB9436WtE7IqIPenlR4BKSaN7s69Z0s4+cRQAL6z2UYEVp1wEwUJgiqTJkqqAa4CHMjtIOkGS0stz0p+7tTf7miXtAw3DqRlU4SCwolWR7RtERLukm4BHgXLgroh4VdIN6e23A58CPi+pHdgPXBMRAXS7b7Y1meVSWZmYM7nOQWBFK+sggHdP9zzSpe32jOXvAt/t7b5m+WbO5DqeWLaZzbvbGDN8UNLlmOWU7yw264U5k1PjBAtXe/pKKz4OArNemDauhiFV5bzgB9BZEXIQmPVCZXkZZ04cyfMeJ7Ai5CAw66UzJ45k+abd7DnQnnQpZjnlIDDrpTMmjCQCXl67I+lSzHLKQWDWS7OaRgDw4hoPGFtxcRCY9VLt4EpOHjOMF9fsSLoUs5xyEJj1wRkTRvDi2h2k7oc0Kw4OArM+OH3CSLbtPcjbW/0AOiseDgKzPjh9wggAFnucwIqIg8CsD6aMGc6w6gqPE1hRcRCY9UF5mZjVVMuLa31EYMXDQWDWR2dMGMmyjbvZd9A3lllxcBCY9dHpE0ZwuCN4Zd3OpEsxywkHgVkfzWgcAcAr6x0EVhxyEgSSLpe0XNJKSTd3s/33JC1Jv56VNCtj21uSXpH0kqSWXNRj1p/qh1czrnYQS3xEYEUi64lpJJUDtwGXkpqDeKGkhyLitYxuq4GPRMR2SVcA84GzM7ZfGBFbsq3FbKDMGF/rIwIrGrk4IpgDrIyIVRFxELgXmJvZISKejYjOyyyeIzVJvVnBmjl+BKu37GXn/kNJl2KWtVwEQSOwNmN9XbqtJ38E/CJjPYDHJC2SNK+nnSTNk9QiqaW1tTWrgs2yNaOxFoBXfVRgRSAXQaBu2rp9EIukC0kFwVcyms+PiDOAK4AbJX24u30jYn5ENEdEc319fbY1m2WlMwiWOAisCOQiCNYBTRnr44ENXTtJmgncAcyNiHfn+4uIDemfm4EHSJ1qMstrI4dW0VQ32JeQWlHIRRAsBKZImiypCrgGeCizg6QJwM+A6yLijYz2oZKGdy4DlwFLc1CTWb+bOX4ES9bvSLoMs6xlfdVQRLRLugl4FCgH7oqIVyXdkN5+O/DXwCjg/0oCaI+IZqABeCDdVgH8S0T8MtuazAbCzMZaHl6yke17DzJyaFXS5Zgdt6yDACAiHgEe6dJ2e8by54DPdbPfKmBW13azQjBjfGqc4JX1O/nwKR63ssLlO4vNjtP0xveCwKyQOQjMjlPNoEpOHD2UJet2JF2KWVYcBGZZmDG+1lcOWcFzEJhlYUZjLRt2ttG6+0DSpZgdNweBWRY6byxb6nECK2AOArMsTGusRfKAsRU2B4FZFoZVVzB59FAHgRU0B4FZlmY01vrUkBU0B4FZlmY01rJxZxtb9njA2AqTg8AsS76xzAqdg8AsS9PG1QCw1PcTWIFyEJhlaXj6DmMfEVihchCY5cB0DxhbAXMQmOVA5x3GWz1gbAXIQWCWAx4wtkLmIDDLgWmN6QFjB4EVoJwEgaTLJS2XtFLSzd1sl6TvpLcvkXRGb/c1KwQ1gyp9h7EVrKyDQFI5cBtwBTAVuFbS1C7drgCmpF/zgO/1YV+zgpAaMN6VdBlmfZaLI4I5wMqIWBURB4F7gbld+swF/jlSngNGSBrby33NCsKMxhrW79jPtr0Hky7FrE9yEQSNwNqM9XXptt706c2+AEiaJ6lFUktra2vWRZvlmgeMrVDlIgjUTVv0sk9v9k01RsyPiOaIaK6v90Thln+me24CK1AVOXiPdUBTxvp4YEMv+1T1Yl+zglAzqJJJo4Z46korOLk4IlgITJE0WVIVcA3wUJc+DwF/kL566BxgZ0Rs7OW+ZgVjemOtTw1Zwck6CCKiHbgJeBRYBtwXEa9KukHSDelujwCrgJXA94E/Ptq+2dZklpQZjbWs37Gf7R4wtgKSi1NDRMQjpL7sM9tuz1gO4Mbe7mtWqGZkDBh/+BSPZVlh8J3FZjk0zVcOWQFyEJjlUO3gSiaOGuIrh6ygOAjMcswDxlZoHARmOTajsZZ12z1gbIXDQWCWYzM8TmAFxkFglmPTxzkIrLA4CMxyrHZIJRPqPGBshcNBYNYPZnjA2AqIg8CsH0z3gLEVEAeBWT/oHDBeusFHBZb/HARm/cBXDlkhcRCY9QMPGFshcRCY9RMPGFuhcBCY9ZPpjbWs3bafHfs8YGz5zUFg1k/eHTBevyvhSsyOzkFg1k+mN9YAHjC2/JdVEEiqk/S4pBXpnyO76dMk6VeSlkl6VdKfZmz7W0nrJb2Ufl2ZTT1m+WTEkCqa6gZ7wNjyXrZHBDcDCyJiCrAgvd5VO/DliDgNOAe4UdLUjO3fjIjZ6ZdnKrOi4gFjKwTZBsFc4O708t3Ax7t2iIiNEbE4vbyb1NzEjVl+rllBmN5Yy5pt+9i571DSpZj1KNsgaIiIjZD6wgfGHK2zpEnA6cDzGc03SVoi6a7uTi1l7DtPUoukltbW1izLNhsYvsPYCsExg0DSE5KWdvOa25cPkjQMuB/4YkR0XkbxPeAkYDawEfiHnvaPiPkR0RwRzfX1nhTcCoMfSW2FoOJYHSLikp62SdokaWxEbJQ0FtjcQ79KUiFwT0T8LOO9N2X0+T7w730p3izfjRxaxfiRgx0ElteyPTX0EHB9evl64MGuHSQJuBNYFhHf6LJtbMbqJ4ClWdZjlndmNNb6yiHLa9kGwa3ApZJWAJem15E0TlLnFUDnA9cBF3VzmejXJL0iaQlwIfClLOsxyzvTG2t5e+s+du73gLHlp2OeGjqaiNgKXNxN+wbgyvTyM4B62P+6bD7frBB0Dhi/un4n5508OuFqzI7kO4vN+pkfSW35zkFg1s9GDq2icYQHjC1/OQjMBoAHjC2fOQjMBsCM8bW8tXUfu9o8YGz5x0FgNgCmd95hvM5HBZZ/HARmA2D2+BEALF6zPdlCzLrhIDAbALVDKjmlYRgL33IQWP5xEJgNkOZJdSxes53DHZF0KWbv4yAwGyDNE0eyu62dNzbtTroUs/dxEJgNkLMm1QHQ8rZPD1l+cRCYDZDxIwczZng1LW9tS7oUs/dxEJgNEEmcNamOFg8YW55xEJgNoOZJI1m/Yz8bduxPuhSzdzkIzAZQ80SPE1j+cRCYDaDTxg5nSFW5xwksr2QVBJLqJD0uaUX6Z7eTz0t6Kz0BzUuSWvq6v1mxqCgv48yJI3l+lYPA8ke2RwQ3AwsiYgqwIL3ekwsjYnZENB/n/mZF4byTRrN8025adx9IuhQzIPsgmAvcnV6+G/j4AO9vVnDOO2kUAL9dtTXhSsxSsg2ChojYCJD+OaaHfgE8JmmRpHnHsT+S5klqkdTS2tqaZdlmyZneWMvwQRX89s0tSZdiBvRizmJJTwAndLPpq334nPMjYoOkMcDjkl6PiKf6sD8RMR+YD9Dc3OyHtVjBKi8T55w4it+s9BGB5YdjBkFEXNLTNkmbJI2NiI2SxgKbe3iPDemfmyU9AMwBngJ6tb9ZsTnvpFE8/tom1m7bR1PdkKTLsRKX7amhh4Dr08vXAw927SBpqKThncvAZcDS3u5vVozOP3k0AL9900cFlrxsg+BW4FJJK4BL0+tIGifpkXSfBuAZSS8DLwAPR8Qvj7a/WbGbMmYYo4dV8azHCSwPHPPU0NFExFbg4m7aNwBXppdXAbP6sr9ZsZPE+SeP5pmVW+joCMrKlHRJVsJ8Z7FZQi78wBi27DnIK+s9j7Ely0FglpAPn1KPBP/xuq+RsGQ5CMwSUje0itObRvDr5Q4CS5aDwCxBF506hpfX7fTjJixRDgKzBF14aupmeh8VWJIcBGYJmjq2hoaaan7lILAEOQjMEiSJi05t4MnlrbQdOpx0OVaiHARmCbtyxgnsPXiYJ9/wwxQtGQ4Cs4Sde+IoRg6p5JFXNiZdipUoB4FZwirKy/idaSewYNlmnx6yRDgIzPLAlTPGsudAO0/59JAlwEFglgfOPWkUI4ZU8rBPD1kCHARmeaCyvIwrZ4zl0VffYXfboaTLsRLjIDDLE586czxthzp4eImPCmxgOQjM8sTpTSM4qX4oP1m0LulSrMQ4CMzyhCQ+3dzEore382brnqTLsRKSVRBIqpP0uKQV6Z8ju+nzAUkvZbx2SfpietvfSlqfse3KbOoxK3T/6fRGysvEfQvXJl2KlZBsjwhuBhZExBRgQXr9fSJieUTMjojZwJnAPuCBjC7f7NweEY903d+slIypGcRlUxu4d+Fa9h/0PQU2MLINgrnA3enlu4GPH6P/xcCbEfF2lp9rVrQ+e/5kdu4/xAMvrk+6FCsR2QZBQ0RsBEj/HHOM/tcAP+7SdpOkJZLu6u7UUidJ8yS1SGppbfVNN1a8zpo0kmnjavjBb1YTEUmXYyXgmEEg6QlJS7t5ze3LB0mqAj4G/CSj+XvAScBsYCPwDz3tHxHzI6I5Iprr6+v78tFmBUUSnz1/Mis27+GpFVuSLsdKwDGDICIuiYjp3bweBDZJGguQ/nm0h6pfASyOiE0Z770pIg5HRAfwfWBOdv85ZsXhqlljGVs7iH9csMJHBdbvsj019BBwfXr5euDBo/S9li6nhTpDJO0TwNIs6zErCtUV5fzxBSfR8vZ2frNya9LlWJHLNghuBS6VtAK4NL2OpHGS3r0CSNKQ9Pafddn/a5JekbQEuBD4Upb1mBWNq89qYmztIL71xBs+KrB+VZHNzhGxldSVQF3bNwBXZqzvA0Z10++6bD7frJhVV5Rz44Un89/+bSmPv7aJy6adkHRJVqR8Z7FZHvvMWU1MGTOM//HwMs9VYP3GQWCWxyrLy/ibq6axZts+7nxmddLlWJFyEJjluQ9OGc3vTGvgH/9jBau37E26HCtCDgKzAnDL3OlUV5Tz5fte4nCHB44ttxwEZgWgoWYQt8ydxuI1O7j9yTeTLseKjIPArEB8bNY4PjpzLP/w2HKe8R3HlkMOArMCIYn/88mZnDxmGDf9eDFrtu5LuiQrEg4CswIytLqC+dc1EwHX3fU8m3e1JV2SFQEHgVmBmTR6KD/47Fm07j7AdXe+wJY9B5IuyQqcg8CsAJ0xYSR3/EEzb2/byye/9yxv+bJSy4KDwKxAnXfyaO753Dns2n+IT37vWQ8g23FzEJgVsDMnjuT+z59H3dAqrrvreb7+6OscaPejKKxvHARmBe7E+mE8eNP5fPrM8dz2qze54ltP8/QKz+JnvecgMCsCQ6oq+NqnZvHDz57F4Qiuu/MFrp3/HM++ucWPsLZjUiH+kjQ3N0dLS0vSZZjlpbZDh/nRc2/zT0+tonX3AU49YTifOauJubMbqRtalXR5liBJiyKi+Yh2B4FZcWo7dJj7F6/j3hfW8sr6nZQJmifWcdFpYzj/pNGcOnY4leU+KVBK+iUIJH0a+FvgNGBORHT77SzpcuDbQDlwR0R0zmRWB/wrMAl4C7g6IrYf63MdBGZ989qGXfxi6UYWLNvMaxt3ATCosoyZjSOY1ljDSfXD0q+hjB5WTVmZEq7Y+kN/BcFpQAfwT8CfdRcEksqBN0hNVbkOWAhcGxGvSfoasC0ibpV0MzAyIr5yrM91EJgdv3d2trHwrW0sXrOdxWt28MY7u9mfMelNZbkYM3wQDTXVnFA7iFFDqxk+qIKawZXUDKqkZnAFwwdVMriynOqKMqory6iuKKeqoiy1XlFGVUUZ5RLlZUJyqOSLnoIg26kql6Xf/Gjd5gArI2JVuu+9wFzgtfTPC9L97gZ+DRwzCMzs+J1QO4irZo3jqlnjAOjoCDbuamNV6x5Wte5l4842Nu9q451dbbz+zm62793K7rZ22o/z8dcSlEuUlYmyjOXyMr23rNS29/Y58julsylzk1A3bUe+h45Y6L5fIfhfn5jBnMl1OX3PrIKglxqBtRnr64Cz08sNEbERICI2ShrT05tImgfMA5gwYUI/lWpWesrKROOIwTSOGMyHptR32yci2H/oMLv2t7Or7RC72w7RdqiDA+2HOXCogwPtHRxsT6+3p9YPdwSHO4KOSP08HEEER7R3/uw8ORHvfmbG59Nl4/v6RTdtvetH4Q2RMrS6POfvecwgkPQE0N2s2V+NiAd78RndxW2f//gjYj4wH1Knhvq6v5kdP0kMqapgSFUFJ9QOSrocy7FjBkFEXJLlZ6wDmjLWxwMb0subJI1NHw2MBTZn+VlmZtZHA3Ht2EJgiqTJkqqAa4CH0tseAq5PL18P9OYIw8zMciirIJD0CUnrgHOBhyU9mm4fJ+kRgIhoB24CHgWWAfdFxKvpt7gVuFTSClJXFd2aTT1mZtZ3vqHMzKxE9HT5qG8rNDMrcQ4CM7MS5yAwMytxDgIzsxJXkIPFklqBt49z99FAPs7p57r6xnX1jevqm3ytC7KrbWJEHHH7eEEGQTYktXQ3ap4019U3rqtvXFff5Gtd0D+1+dSQmVmJcxCYmZW4UgyC+UkX0APX1Teuq29cV9/ka13QD7WV3BiBmZm9XykeEZiZWQYHgZlZiSvJIJA0W9Jzkl6S1CJpTtI1dZL0BUnLJb2antM5b0j6M0khaXTStQBI+rqk1yUtkfSApBEJ13N5+v/dyvQc3ImT1CTpV5KWpX+n/jTpmjJJKpf0oqR/T7qWTpJGSPpp+ndrmaRzk64JQNKX0v8Pl0r6saSczRBUkkEAfA347xExG/jr9HriJF1Iah7nmRExDfj7hEt6l6QmUo8KX5N0LRkeB6ZHxEzgDeAvkipEUjlwG3AFMBW4VtLUpOrJ0A58OSJOA84BbsyTujr9KanH0+eTbwO/jIhTgVnkQX2SGoE/AZojYjpQTmpul5wo1SAIoCa9XMt7M6Yl7fPArRFxACAi8mnGtm8Cf04ezfIaEY+l57sAeI7U7HdJmQOsjIhVEXEQuJdUqCcqIjZGxOL08m5SX2qNyVaVImk88LvAHUnX0klSDfBh4E6AiDgYETsSLeo9FcBgSRXAEHL4vVWqQfBF4OuS1pL6V3di/5Ls4hTgQ5Kel/SkpLOSLghA0seA9RHxctK1HMUfAr9I8PMbgbUZ6+vIky/cTpImAacDzydcSqdvkfrHRUfCdWQ6EWgFfpA+ZXWHpKFJFxUR60l9V60BNgI7I+KxXL3/MecsLlSSngBO6GbTV4GLgS9FxP2SriaV/tnOzZyLuiqAkaQO4c8C7pN0YgzANb7HqOsvgcv6u4buHK2uiHgw3eerpE6B3DOQtXWhbtry5uhJ0jDgfuCLEbErD+r5KLA5IhZJuiDhcjJVAGcAX4iI5yV9G7gZ+Kski5I0ktQR5mRgB/ATSb8fET/KxfsXbRBERI9f7JL+mdS5SYCfMICHpseo6/PAz9Jf/C9I6iD1gKnWpOqSNIPUL9/LkiB1+mWxpDkR8U5SdWXUdz3wUeDigQjMo1gHNGWsjydPTjlKqiQVAvdExM+SriftfOBjkq4EBgE1kn4UEb+fcF3rgHUR0XnU9FNSQZC0S4DVEdEKIOlnwHlAToKgVE8NbQA+kl6+CFiRYC2Z/o1UPUg6Bagi4ScgRsQrETEmIiZFxCRSf1HOGIgQOBZJlwNfAT4WEfsSLmchMEXSZElVpAbyHkq4JpRK7zuBZRHxjaTr6RQRfxER49O/U9cA/5EHIUD693qtpA+kmy4GXkuwpE5rgHMkDUn/P72YHA5iF+0RwTH8F+Db6UGXNmBewvV0ugu4S9JS4CBwfcL/ys133wWqgcfTRyvPRcQNSRQSEe2SbgIeJXVFx10R8WoStXRxPnAd8Iqkl9JtfxkRjyRXUt77AnBPOtBXAZ9NuB7Sp6l+CiwmdRr0RXL4qAk/YsLMrMSV6qkhMzNLcxCYmZU4B4GZWYlzEJiZlTgHgZlZiXMQmJmVOAeBmVmJ+/+hdzkxbUCUdwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "inps = np.linspace(-8, 8, 10000)\n",
    "plt.plot(inps, -np.tanh(inps))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "7ddff01a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# integral of tanh is ln |cosh (x)|\n",
    "def make_g(f, transition_pt):\n",
    "    return lambda J, x: -np.log(np.abs(np.cosh(f(J, x) - transition_pt)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "633811e6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.9464173888754903\n",
      "-0.9464173888754903 -0.8545211150237193\n",
      "-1.00693663436774 -0.9303740341919693\n",
      "-1.00693663436774 -0.9290480898417486\n",
      "-1.00693663436774 -0.9350048791113655\n",
      "-1.00693663436774 -0.9057240012752328\n",
      "-1.00693663436774 -0.9290480898417486\n",
      "-1.00693663436774 -0.9350048791113655\n",
      "-1.00693663436774 -0.9057240012752328\n",
      "-1.00693663436774 -0.9290480898417486\n",
      "-1.00693663436774 -0.9350048791113655\n"
     ]
    }
   ],
   "source": [
    "N = 100\n",
    "J = np.random.normal(size=(N,N))\n",
    "\n",
    "x = random_state()\n",
    "x, best = local_search(J, cost, x)\n",
    "x_best = x\n",
    "print(best)\n",
    "for i in range(100):\n",
    "    g = make_g(cost, best-0.01) # the 0.01 controls the offset\n",
    "    x, _ = local_search(J, g, x)\n",
    "    x, test = local_search(J, cost, x)\n",
    "    if i % 10 == 0:\n",
    "        print(best, test)\n",
    "    if test < best:\n",
    "        best = test\n",
    "        x_best = x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "70ac1351",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1., -1.,  1., -1., -1.,  1., -1.,  1., -1., -1.,  1., -1.,  1.,\n",
       "        1., -1., -1.,  1., -1.,  1.,  1.,  1., -1., -1., -1., -1., -1.,\n",
       "        1.,  1., -1., -1.,  1.,  1., -1., -1., -1.,  1., -1.,  1., -1.,\n",
       "       -1., -1.,  1.,  1., -1., -1.,  1., -1.,  1., -1., -1., -1.,  1.,\n",
       "        1., -1.,  1.,  1., -1.,  1.,  1., -1.,  1.,  1.,  1.,  1., -1.,\n",
       "        1.,  1., -1.,  1., -1., -1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,\n",
       "        1.,  1.,  1., -1., -1., -1.,  1.,  1., -1., -1.,  1., -1.,  1.,\n",
       "        1.,  1.,  1.,  1., -1., -1.,  1.,  1., -1.])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_best"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f473eb9",
   "metadata": {},
   "source": [
    "Other things we could try:\n",
    "* change the width of the tanh function\n",
    "* implement the other functions like hBOA or GA ? (These are quite complicated...)\n",
    "* test at small N? (implement branch n bround?) or not, it's cool... might want to see if we're getting the right minimum...."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ebfb994",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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
}