guanjihuan/GitHub_collection_pykan
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
GitHub_collection_pykan/tutorials/Example_6_PDE_interpretation.ipynb
Ziming Liu 830418878c
updated tutorials
2024-08-11 13:02:16 -04:00

321 lines
48 KiB
Plaintext
Raw Blame History

{
"cells": [
{
"cell_type": "markdown",
"id": "5d904dee",
"metadata": {},
"source": [
"# Example 6: Solving Partial Differential Equation (PDE)"
]
},
{
"cell_type": "markdown",
"id": "7d568912",
"metadata": {},
"source": [
"We aim to solve a 2D poisson equation $\\nabla^2 f(x,y) = -2\\pi^2{\\rm sin}(\\pi x){\\rm sin}(\\pi y)$, with boundary condition $f(-1,y)=f(1,y)=f(x,-1)=f(x,1)=0$. The ground truth solution is $f(x,y)={\\rm sin}(\\pi x){\\rm sin}(\\pi y)$."
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "0e2bc449",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"checkpoint directory created: ./model\n",
"saving model version 0.0\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"pde loss: 2.23e+00 | bc loss: 9.47e-03 | l2: 3.63e-03 : 100%|███████| 20/20 [00:24<00:00, 1.25s/it]\n"
]
}
],
"source": [
"from kan import *\n",
"import matplotlib.pyplot as plt\n",
"from torch import autograd\n",
"from tqdm import tqdm\n",
"\n",
"dim = 2\n",
"np_i = 21 # number of interior points (along each dimension)\n",
"np_b = 21 # number of boundary points (along each dimension)\n",
"ranges = [-1, 1]\n",
"\n",
"model = KAN(width=[2,2,1], grid=5, k=3, seed=4)\n",
"\n",
"def batch_jacobian(func, x, create_graph=False):\n",
" # x in shape (Batch, Length)\n",
" def _func_sum(x):\n",
" return func(x).sum(dim=0)\n",
" return autograd.functional.jacobian(_func_sum, x, create_graph=create_graph).permute(1,0,2)\n",
"\n",
"# define solution\n",
"sol_fun = lambda x: torch.sin(torch.pi*x[:,[0]])*torch.sin(torch.pi*x[:,[1]])\n",
"source_fun = lambda x: -2*torch.pi**2 * torch.sin(torch.pi*x[:,[0]])*torch.sin(torch.pi*x[:,[1]])\n",
"\n",
"# interior\n",
"sampling_mode = 'random' # 'radnom' or 'mesh'\n",
"\n",
"x_mesh = torch.linspace(ranges[0],ranges[1],steps=np_i)\n",
"y_mesh = torch.linspace(ranges[0],ranges[1],steps=np_i)\n",
"X, Y = torch.meshgrid(x_mesh, y_mesh, indexing=\"ij\")\n",
"if sampling_mode == 'mesh':\n",
" #mesh\n",
" x_i = torch.stack([X.reshape(-1,), Y.reshape(-1,)]).permute(1,0)\n",
"else:\n",
" #random\n",
" x_i = torch.rand((np_i**2,2))*2-1\n",
"\n",
"# boundary, 4 sides\n",
"helper = lambda X, Y: torch.stack([X.reshape(-1,), Y.reshape(-1,)]).permute(1,0)\n",
"xb1 = helper(X[0], Y[0])\n",
"xb2 = helper(X[-1], Y[0])\n",
"xb3 = helper(X[:,0], Y[:,0])\n",
"xb4 = helper(X[:,0], Y[:,-1])\n",
"x_b = torch.cat([xb1, xb2, xb3, xb4], dim=0)\n",
"\n",
"steps = 20\n",
"alpha = 0.01\n",
"log = 1\n",
"\n",
"def train():\n",
" optimizer = LBFGS(model.parameters(), lr=1, history_size=10, line_search_fn=\"strong_wolfe\", tolerance_grad=1e-32, tolerance_change=1e-32, tolerance_ys=1e-32)\n",
"\n",
" pbar = tqdm(range(steps), desc='description', ncols=100)\n",
"\n",
" for _ in pbar:\n",
" def closure():\n",
" global pde_loss, bc_loss\n",
" optimizer.zero_grad()\n",
" # interior loss\n",
" sol = sol_fun(x_i)\n",
" sol_D1_fun = lambda x: batch_jacobian(model, x, create_graph=True)[:,0,:]\n",
" sol_D1 = sol_D1_fun(x_i)\n",
" sol_D2 = batch_jacobian(sol_D1_fun, x_i, create_graph=True)[:,:,:]\n",
" lap = torch.sum(torch.diagonal(sol_D2, dim1=1, dim2=2), dim=1, keepdim=True)\n",
" source = source_fun(x_i)\n",
" pde_loss = torch.mean((lap - source)**2)\n",
"\n",
" # boundary loss\n",
" bc_true = sol_fun(x_b)\n",
" bc_pred = model(x_b)\n",
" bc_loss = torch.mean((bc_pred-bc_true)**2)\n",
"\n",
" loss = alpha * pde_loss + bc_loss\n",
" loss.backward()\n",
" return loss\n",
"\n",
" if _ % 5 == 0 and _ < 50:\n",
" model.update_grid_from_samples(x_i)\n",
"\n",
" optimizer.step(closure)\n",
" sol = sol_fun(x_i)\n",
" loss = alpha * pde_loss + bc_loss\n",
" l2 = torch.mean((model(x_i) - sol)**2)\n",
"\n",
" if _ % log == 0:\n",
" pbar.set_description(\"pde loss: %.2e | bc loss: %.2e | l2: %.2e \" % (pde_loss.cpu().detach().numpy(), bc_loss.cpu().detach().numpy(), l2.detach().numpy()))\n",
"\n",
"train()"
]
},
{
"cell_type": "markdown",
"id": "e2246bab",
"metadata": {},
"source": [
"Plot the trained KAN"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "02e2a0ba",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAFICAYAAACcDrP3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABO8klEQVR4nO3dd1hT59sH8O9Jwh6yQQSUJQi4657UWXfVOqvVqnVVq9XWOuqedVFHa2tr66iram21irYWF+6JggKCoCyRESQkQJLzvH/4S15wIwkngftzXVxXawY3hPt8zznPc57DMcYYCCGEEB0SCV0AIYSQyofChRBCiM5RuBBCCNE5ChdCCCE6R+FCCCFE5yhcCCGE6ByFCyGEEJ2jcCGEEKJzFC6EEEJ0jsKFEEKIzlG4EEII0TkKF0IIITpH4UIIIUTnKFwIIYToHIULIYQQnZMIXQAhxoAxhuzsbMhkMlhbW8PR0REcxwldFiEGi45cCHkFqVSKb7/9Fv7+/nB2doa3tzecnZ3h7++Pb7/9FlKpVOgSCTFIHN2JkpAXO3bsGPr16we5XA7g6dGLhuaoxdLSEvv370eXLl0EqZEQQ0XhQsgLHDt2DN27dwdjDDzPv/R5IpEIHMfh77//poAhpAQKF0KeIZVK4eHhAYVC8cpg0RCJRLCwsEBKSgrs7Oz0XyAhRoDGXAh5xtatWyGXy98oWACA53nI5XJs27ZNz5URYjzoyIWQEhhj8Pf3R2JiIsrSGhzHwcfHB/Hx8TSLjBBQuBBSSlZWFpydncv1ekdHRx1WRIhxotNihJQgk8nK9fr8/HwdVUKIcaNwIaQEa2vrcr3exsZGR5UQYtwoXAgpwdHREb6+vmUeN+E4Dr6+vnBwcNBTZYQYFwoXQkrgOA6TJk16q9dOnjyZBvMJ+R8a0CfkGXSdCyHlR0cuhDzDzs4O+/fvB8dxEIle3SKaK/QPHDhAwUJICRQuhLxAly5d8Pfff8PCwgIcxz13ukvzbxYWFjhy5Ag6d+4sUKWEGCYKF0JeokuXLkhJSUFYWBh8fHxKPebj44OwsDCkpqZSsBDyAjTmQsgbYIwhIiICHTp0wIkTJxAaGkqD94S8Ah25EPIGOI7TjqnY2dlRsBDyGhQuhBBCdI7ChRBCiM5RuBBCCNE5ChdCCCE6R+FCCCFE5yhcCCGE6ByFCyGEEJ2jcCGEEKJzFC6EEEJ0jsKFEEKIzlG4EEII0TkKF0IIITpH4UIIIUTnKFwIIYToHIULIYQQnaNwIYQQonMULoS8hlKpRGpqKu7cuQMASEhIQE5ODnieF7gyQgwX3eaYkJeQSqXYv38/fvvtN0RHRyM/Px/FxcUwNzeHs7Mz2rRpg1GjRqFVq1aQSCRCl0uIQaFwIeQFzp8/j6lTpyIqKgpNmjRB9+7dUa9ePVhbW0MqleLq1as4dOgQ7t27h4EDB2Lx4sVwdnYWumxCDAaFCyHPOH78OEaMGAFra2ssW7YM3bp1Q3FxMXbv3o2ioiLY2tpi0KBBUCqV2L17N+bPn4/g4GBs374drq6uQpdPiEGgcCGkhLi4OHTt2hVWVlbYvXs3goKCwHEcEhMT0ahRI+Tl5cHb2xtXr16Fvb09GGM4e/YshgwZgvbt2+Onn36CmZmZ0D8GIYKjAX1C/ketVmPp0qXIzc3Fhg0btMHyKhzHoXXr1vjmm2/w559/Ijw8vIKqJcSwUbgQ8j/37t3DoUOH0LdvX7Ru3fq1waLBcRz69OmD5s2bY/PmzVCpVHqulBDDR1NcCPmfc+fOQSaToV+/fkhKSkJBQYH2sZSUFKjVagBAcXExoqOjYWtrq33c3d0dffv2xfz585GRkQEPD48Kr58QQ0LhQsj/3L17F5aWlvDx8cHYsWMRGRmpfYwxhqKiIgBAWloaOnXqpH2M4zisXr0adevWhVwuR1paGoULqfIoXAj5H4VCAYlEAjMzMxQVFaGwsPCFz2OMPfeYSqWChYVFqRAipCqjcCFV3v379xEREYHTp09DLpdDKpWiWbNmsLKy0j5HoVDg3Llz2hBp2bKl9sJJjuPg5eWFzMxMqFQqxMfHo0mTJjA3NxfqRyJEcDQVmVQ5Dx48wMmTJxEREYGIiAgkJyeD4zh4e3sjOTkZGzduxOjRo0u9JjExEU2aNEFeXh5q1aqFK1euwM7OTvs4x3GYNWsWVq5cCZ7nYWZmhhYtWqB9+/YIDQ1Fs2bNaIoyqVIoXEill5qaioiICG2gJCYmAgDq16+v3fi3bdsWPM+jdevWsLe3R3h4eKkB+5dd5wI8PU2WlpaGdu3aoWfPnvjoo49w8uRJnDx5EqdOnYJUKtUe7Wi+X5MmTWBqairI74OQikDhQiqdjIyMUmESHx8PAAgJCdFu3Nu1awdHR8fnXrtx40ZMmzYNc+bMwVdffaU99fWqcCksLMSUKVNw6NAh/PfffwgICNC+n1qtxs2bN7W1nD59Gk+ePIGlpSVatWqF0NBQhIaGonHjxjAxMamA3w4hFYPChRi9zMxM7ZFCREQE7t69CwCoU6dOqTBxcXF57XsVFBTg448/xpEjR7BgwQKMHz8e5ubmuH//Ppo2bao9LXbp0iXY2dkhPz8fS5YswQ8//IC1a9di5MiRr3x/lUqF69eva8PvzJkzkMlksLa2RuvWrbVh07BhQ1oMkxg1ChdidLKysnDq1CltmERHRwMAateurQ2T9u3bw83N7a3e//Hjx5g4cSIOHz6MLl26YOrUqahTpw5iY2PB8zxMTU3h5+eHS5cuYdWqVbhx4wYWLlyI8ePHQywWl+l7KZVKXL16VRs2Z8+ehVwuh62tLdq0aaMNm/r165f5vQkREoULMXi5ubk4deqUdgMcFRUFAPD19S0VJjVq1NDZ9ywoKMDmzZuxbt06PHr0CD4+PvD394eNjQ1yc3MRGxuLtLQ0NG7cGPPmzUO7du0gEpV/wYvi4mJcvnxZG5yRkZEoLCyEnZ0d2rZtqw2bunXr6uT7EaIvFC7E4OTl5eH06dPaMLlx4wYYY6hVq5Y2SEJDQ+Hp6an3WjIyMnDixAmcOnUKiYmJKCwshL29PUJCQtC5c2c0a9YMlpaWevv+RUVFuHjxojZszp8/j6KiIjg4OKBdu3basAkODn7j5WoIqQgULkRw+fn5OHPmjHYDeu3aNfA8Dw8PD+3GMzQ0FLVq1RK0TrVaDcYYRCKRYEcNhYWFOH/+vPZ3deHCBSiVSjg7O5cKm8DAQAobIigKF1LhZDIZIiMjtUcmV65cgVqthru7e6kjEx8fH9pAvoZcLse5c+e0YXPp0iWoVCq4urpqf4+hoaHw9/en3yWpUBQuRO80G0BNmJTcAJYME9oAlp9MJtP+riMiIkoFd8mwoeAm+kbhQnROc+pGEyYlT920b99eu5GjUzf69+TJE+1RYslTjp6enqXCRuhTjqTyoXAh5aYZdNaESclB55JhQoPOwpNKpaXGt0pOligZNhUxWYJUbhQupMw002U1e8Pnzp3TTpdt166ddiNF02UNX05ODs6cOaP9LDXTvH18fEpNpnB3dxe4UmJsKFzIa5W80E9z7YXmQj/NtRft27enC/0qgaysLO008JIXqPr7+2uDpjwXqJKqg8KFPKfkEiURERE4e/asdomSNm3aaI9MaImSyi8zM1N7AWvJpXUCAwNLhY2zs7PAlRJDQ+FCtIsrajYgZ86c0S6uqFnvqn379rS4IkF6enqpsNEsChocHKwNm5ctCkqqFgqXKojnedy6dUu7gTh9+jSkUinMzc21K/W2b9+eloUnr5Wamlrq3jia2xnUq1dPGzZt27bVriBNqg4KlyqAMYbo6GjtBuDUqVPIycnR3tBKEyZ0QytSXi+7EVuDBg20YdOmTRtUq1ZN6FKJnlG4VEKMMdy9e7dUmDx+/BgmJiZo3ry5NkxatGhBt+IlenX//v1SYZOSkgKRSIRGjRppw6Z169awsbERulSiYxQulQBjDPHx8doGPnnyJB49egSJRIKmTZtqw6Rly5Z6XWSRkFdhjCExMVH7dxoREYH09HSIxWK888472rBp1aoVrKyshC6XlBOFixF6tklPnjyJtLS055q0ZcuWsLa2FrpcQl7oTXaKQkND0aJFC9opMkIULkYiKSmpVBM+fPjwudMLrVq1KnXfd0KMybOnc0+ePImsrCyYmpqiWbNmdDrXyFC4GInAwEDExcWVGhht3bo17OzshC6NEL3geR4xMTHPTUTZvXs3Bg4cKHR55DUoXAghhOgcLfxECCFE52jtDh3RDE5mZ2cLXUq5iEQihISE0GwdUmbUA6QkOi2mIzzPY8KECfDw8ICFhQXUajVMTU2Nbon506dPY968eahXr57QpRAjo+kBT09Po5ulqFQqoVAoYGNjgzNnzlAP6AAdueiQmZkZmjdvjqVLl8LExAQ7d+40qjWWGGOQyWSg/Q3ytszMzDB69Gi4uroKXcobY4zh2LFjmD59Ovr164fg4GDqAR2gMRcdUyqVOH/+PK5cuQKZTCZ0OYSQ1+B5Hrt27UJ0dDRu3rxJt43QEQoXHbO2toapqSmKioooXAgxAmlpafj3338hFosxcOBAuo2EjlC46JiVlRXMzc2hVCqRn58vdDmEkFdgjCE8PBzp6enw9PREaGio0CVVGhQuOmZhYaEd0M/LyxO6HELIKxQVFWHPnj1gjKF79+5wcXERuqRKg8JFx0xNTWFtbQ2e5yGVSoUuhxDyEowx3L59GxcvXoS5uTkGDBhgdLM7DRmFi46ZmJjAxsYGjDHk5uYKXQ4h5BV+//13yGQy1K9fH40bNxa6nEqFwkXHJBKJdr2v7OxsmtJIiIHKysrCH3/8AY7jMGjQIFp5WcdoWoSOiUQiODg4AECFXKnMGENRURHi4uJw48YNPH78GPb29ggMDIS/vz8cHR3BcRwd7hNSAmMMJ06cQGJiIlxcXNCrVy9wHEc7gzpE4aIHFREujDGoVCr8888/WL16NS5dulRq6rOpqSnc3d3Rt29fjB07Fn5+fhCJ6ECVEODp9Wjbtm2DWq1Gly5dULNmTaFLqnQoXHSM4zjtVfm5ubngeV7nF2UxxpCXl4f58+dj8+bNkMvlsLS0RL169eDm5oacnBwkJyfjwYMHWLNmDXbt2oUxY8agTZs2qFWrFpydnWFlZaWti45qSFWiGcg/e/YszMzMMGzYMNrx0gMKFz0oGS5qtVqn4cIYQ3Z2NsaNG4c//vgDYrEYffv2xfTp01GvXj3tNTaZmZk4evQoNm7ciNu3b2PhwoUQi8WwtraGi4sLvL29ERgYiODgYAQEBKBmzZpwcnLSWZ2EGCrGGHbs2IH8/Hw0b94czZs3px0sPaBw0QMnJydwHAepVAqlUglTU1OdvXd+fj4+/fRTHDhwAFZWVli0aBHGjh0Lc3NzbYOIxWJ4eXnhk08+Qb9+/fDDDz8gPDwcDx48QHZ2NhISEhAfH4/jx48DeDoJwcbGBh07dkRQUJDOaiXEEKWnp+PAgQPgOA7Dhg2j1Y/1hMJFDxwcHCASiZCfn4/CwkKd/fEqlUosWbIEv//+OywtLbF27VqMHDnypUdGHMfByckJs2bNwhdffIEnT54gIyMDCQkJiImJQUxMDOLj4/Hw4UPk5ubC1NSUlr4glRpjDAcPHsSDBw/g5eWF3r1701GLntCWRA/s7e0hkUhQUFAAuVyuk5WRGWP4888/sX79eohEIsyaNQsjRox4o1NuHMfB1NQUTk5OcHJyQkhICHr16gXGGJRKJaRSKdLS0mBubo4DBw6Uu1ZCDFV+fj5+/fVXMMbQv39/VK9eXeiSKi0KFz2ws7ODmZkZFAoFnjx5Uu73Y4whOTkZM2fOhEKhwMCBAzFlypRyHWVopiebmZnB1dUVrq6uNA2TVGqa6cc3b96EnZ0dhg8fTkctekRTJPTA2toaFhYWKC4u1skSMCqVCkuXLsW9e/fg6+uLpUuXwsLCovyFElKFFBUV4YcffoBSqUTXrl0RFBRE4aJHFC56YGlpCRsbG6jV6nJf68IYw5kzZ7Bz506YmJhg7ty58Pb2pqYgpAwYYzh//jxOnz4NCwsLfPLJJ3TfFj2jcNEDc3Nz2NnZged5PH78uFzvJZfLsXTpUhQUFKBjx47o378/BQshZaRSqfDdd99BoVCgXbt2aNGiBfWRnlG46IFEItEO4j969OitxzIYYzhy5AhOnz4NGxsbzJw5k06HEVJGjDFcuXIF4eHhMDU1xYQJE2BmZiZ0WZUehYseiMViODs7AwAyMjLe+n0KCgoQFhYGpVKJPn360MVehLwFlUqFdevWQSaToWXLlujQoQP1UQWgcNETzRTHjIyMtzpyYYzh2LFjuHz5MqpVq4bJkyfTNSiElBFjDBcvXsThw4dhamqKKVOm0NF/BaFw0RM3NzcAQGZmJtRqdZlfX1hYiO+//x5KpRI9evRAgwYNaG+LkDIqKirCqlWrIJPJ0LZtW3Tq1In6qIJQuOgBx3GoXr06OI5DdnY2ioqKyvR6zcyWyMhIWFlZYdy4cTSzhZAyYozh33//xbFjx2BhYYHp06fTUUsFonDREzc3N4jFYuTm5qKgoKBMr1Wr1fj5559RWFiItm3bokmTJrS3RUgZPXnyBMuWLUNhYSG6deuG9u3bUx9VIAoXPXF2doaZmRny8/ORl5dXptfGxcUhPDwcEokEo0eP1unCl4RUBYwxbNu2DRcvXoS9vT2++uor6qMKRuGiJ/b29rCyskJhYSGysrLe+HWMMezatQs5OTkIDg7Gu+++S3tbhJQBYwzx8fFYuXIleJ7H6NGj0bBhQ+qjCkbhoic2Njaws7ODUqlEenr6G78uOzsbe/fuBcdxGDp0KKpVq6bHKgmpfIqLi7Fw4UI8fPgQderUwZQpU2jMUgAULnpibm6uXQzy4cOHbzQdmTGGf/75BwkJCXBxcUHfvn1pb4uQMmCM4ffff8e+fftgbm6O+fPn08rHAqFw0ROJRIIaNWoAAJKSkt7oNcXFxdr7er/33nuoVauW/gokpJJhjOHu3buYM2cOioqKMHjwYLpfi4AoXPSE4zjUrFkTAJCcnAye51/5fMYYoqKicPbsWZibm9N9vQkpA8YYnjx5gs8//xzJyckICgrC/PnzYWJiInRpVRZtvfTIx8cHAJCSkoLi4uJXPpcxhp07d0Imk6FRo0Zo1qwZ7XER8oaUSiUWL16M48ePo1q1ali9ejU8PT2phwRE4aInHMehVq1aEIvFyMjIgEwme+XzMzIy8Mcff2jv621paVlBlRJi3DTXhW3YsEF7l1a6El94FC565OHhAQsLC0ilUmRmZr70eYwxHD58GA8fPoSHhwd69uxJjUHIG+B5HgcOHMDMmTNRVFSEjz76CJMmTaLZYQaAwkWPXFxc4OjoCIVCgeTk5Jc+T6FQYNu2beB5Hr1796bZLYS8AZ7ncejQIUyYMAF5eXl47733sGLFCpibmwtdGgGFi17Z2NjAw8MDarUasbGxL5yOrFlH7MqVK7C2tsawYcPoqIWQ11Cr1di9ezdGjRqFrKwstG/fHj/++CMcHByofwwEhYsemZqaIiAgAABw+/btFz5HrVZjy5YtKCoqQtu2bVG/fn1qDkJegjEGhUKB1atXY+zYscjOzsa7776LrVu3wt3dnXrHgNANQvSsbt26AICYmBgUFRWVOmRnjCEmJgZHjx6FiYkJRo0aResfEfISjDE8ePAAs2bNwt69e8HzPPr27YuNGzfC1dWVgsXA0JGLHnEch3r16sHExAT3799/bo0xxhi2bNmC3Nxc1K9fHx07dqQGIeQZjDEUFBRg+/bt6NSpE3bu3AkTExNMmzYNW7ZsoWAxUHTkomf+/v5wcHBATk4OYmNj4eHhAeBpwyQkJGDPnj0QiUQYM2YMbGxsBK6WEMPBGEN+fj5OnDiBdevWITIyEkqlEr6+vli6dCnef/99ukjSgFG46JmLiwsCAwNx6tQpnD17VrvKMc/z2LRpEzIyMhAcHEzriJEqTTPZhed5PHnyBLGxsQgPD8fBgwcRExMDpVIJW1tbjB49GjNmzICXlxf1i4GjcNEzU1NTtG3bFqdOncKJEycwY8YMmJmZITo6Gtu2bYNIJMKECRPg6OgodKmEVAhNkKhUKuTl5SE5ORm3bt1CVFQUYmNjERsbi9TUVBQWFgIAHBwc0L17d3z66ado1KgRxGIxBYsRoHCpAJ07d8bKlStx48YNxMTEICgoCEuWLEFWVhbeeecdDBo0iJqFVFqMMfA8j/z8fKSkpCAmJgZXrlzB9evXERcXh8zMTG2QAIBIJIKtrS0aN26Mnj17onfv3vDz86NQMTIULnrGcRwaNGiA+vXr4+LFi1i/fj3q1auHgwcPwtzcHLNnz4a9vb3QZRKiMzzPQyaTIT09HXfu3MG1a9dw9epV3L17F+np6ZDL5dqjF7FYDFtbWwQEBCAoKAiBgYGoV68eQkJC4OHhATMzMwoUI0XhUgGsrKwwceJEXLt2Ddu3b4dIJIJKpcK4cePQrVs3ah5SqVy4cAHz5s3DgwcPkJ+fr10RXCQSwcrKShskjRs3RqNGjVC7dm24ubnB3NwcHMdRP1QSFC4VgOM49O/fH9euXcNPP/0EnucxePBgLF68mGa7kErH1tYWSUlJUKvVqFmzJgIDA9GoUSM0atQIQUFBcHd3h5WVFUQiEQVJJUbhokOMMeTm5r40MGbMmIHevXuD53nUqVMHAJCbm1uRJb5WyXPfhJQVYwwODg5Yv349atasCU9PT9ja2pZaSFIzkG+oqAd0g8JFRzQ3B1u3bt0brcj6xx9/VEBVZadQKFCtWjWhyyBGSNMDmzdvhlgsxpUrV4Qu6a1QD+gGx97k5u7ktRhjL1yY0hjReW/yNqgHSEkULoQQQnSO1hYjhBCicxQuhBBCdI7CxUhcu3YNHMfh2rVrQpdCiGCoD4wHhQshhBCdo3AhhBCicxQuhBBCdI7ChRBCiM5RuBBCCNE5ChdCCCE6R+FCCCFE5yhcCCGE6ByFCyGEEJ2jcCGEEKJzFC6EEEJ0jsKFEEKIzlG4EEII0TkKF0IIITpH4UIIIUTnKFyMAGMMubm5AIDc3NxKc59yQsqC+sC4ULgYMKlUim+//Rb+/v7o2LEjAKBjx47w9/fHt99+C6lUKmyBhFQA6gPjxDGKf4N07Ngx9OvXD3K5HABK7aVxHAcAsLS0xP79+9GlSxdBaiRE36gPjBeFiwE6duwYunfvDsYYeJ5/6fNEIhE4jsPff/9NjUUqHeoD40bhYmCkUik8PDygUChe2VAaIpEIFhYWSElJgZ2dnf4LJKQCUB8YPxpzMTBbt26FXC5/o4YCAJ7nIZfLsW3bNj1XRkjFoT4wfnTkYkAYY/D390diYmKZZsJwHAcfHx/Ex8drz0MTYqyoDyoHChcDkpWVBWdn53K93tHRUYcVEVLxqA8qBzotZkBkMlm5Xp+fn6+jSggRDvVB5UDhYkCsra3L9XobGxsdVUKIcKgPKgcKFwPi6OgIX1/fMp8v5jgOvr6+cHBw0FNlhFQc6oPKgcLFgHAch0mTJr3VaydPnkyDmKRSoD6oHGhA38DQ/H5CqA8qAzpyMTB2dnbYv38/OI6DSPTqj0dzZfKBAweooUilQn1g/ChcDFCXLl3w999/w8LCAhzHPXeYr/k3CwsLHDlyBJ07dxaoUkL0h/rAuFG4GKguXbogJSUFYWFh8PHxKfWYj48PwsLCkJqaSg1FKjXqA+NFYy5GgDGGiIgIdOjQASdOnEBoaCgNWpIqh/rAuNCRixHgOE57LtnOzo4ailRJ1AfGhcKFEEKIzlG4EEII0TkKF0IIITpH4UIIIUTnKFwIIYToHIULIYQQnaNwIYQQonMULoQQQnSOwoUQQojOUbgQQgjROQoXQgghOkfhQgghROcoXAghhOgchQshhBCdo3AhhBCicxQuhBBCdI7CxcDJZDLExcXh1q1bAICMjAwUFxcLXBUhFUupVCI1NRV37twBACQkJCAnJwc8zwtcGXkZus2xgUpMTMRPP/2Ev/76Cw8fPoRSqURRURFsbW3RsGFDfPTRR+jbty9sbGyELpUQvZFKpdi/fz9+++03REdHIz8/H8XFxTA3N4ezszPatGmDUaNGoVWrVpBIJEKXS0qgcDEwarUau3btwqxZs6BQKPDee++hU6dO8PLyAs/zuHfvHo4ePYqIiAg0atQI69evR1BQkNBlE6Jz58+fx9SpUxEVFYUmTZqge/fuqFevHqytrSGVSnH16lUcOnQI9+7dw8CBA7F48WI4OzsLXTb5HwoXA8LzPDZt2oQvv/wSbdu2xfLlyxEcHIxLly7h8uXLAIDOnTvD19cX58+fx7Rp05Cfn499+/YhJCRE4OoJ0Z3jx49jxIgRsLa2xrJly9CtWzcUFxdj9+7d2iP4QYMGQalUYvfu3Zg/fz6Cg4Oxfft2uLq6Cl0+AQBGDEZERASzs7Nj/fv3Zzk5OYznecYYY3PmzGEAGAC2fft2xhhjPM+z5ORk1rJlS9a6dWuWm5srYOWE6E5sbCzz9vZmISEh7Pbt29o+SEhIYNWqVWMAmLe3N8vJyWGMPe2F06dPMw8PD/bhhx+ywsJCIcsn/0MD+gZCoVBg4cKFcHV1xdq1a2FnZweO4176fI7j4OnpifXr1yMuLg47duyowGoJ0Q+1Wo2lS5ciNzcXGzZsQFBQ0Cv7AHjaC61bt8Y333yDP//8E+Hh4RVULXkVChcDcfXqVVy4cAETJkxAjRo1XttQwNOmatCgAQYMGIBff/0Vcrm8AiolRH/u3buHQ4cOoW/fvmjduvUb9QHwtBf69OmD5s2bY/PmzVCpVHqulLwOTa8wECdPnoSZmRk6duyIO3fulGqOR48eaf/7wYMHiIqK0v6/nZ0d+vTpgx07diApKYkG94lRO3fuHGQyGfr164ekpCQUFBRoH0tJSYFarQYAFBcXIzo6Gra2ttrH3d3d0bdvX8yfPx8ZGRnw8PCo8PrJ/6NwMRCxsbFwcXGBiYkJOnbsiMzMTO1jJYNm/vz5WLRokfb/Bw8ejLlz50IikeDBgwcULsSo3b17F5aWlvDx8cHYsWMRGRmpfYwxhqKiIgBAWloaOnXqpH2M4zisXr0adevWhVwuR1paGoWLwChcDABjDIWFhTAzM4NYLEZhYSEKCwtf+FylUgmlUqn9/+LiYpiamkIkEuHy5cto1KgRXFxcKqp0QnSiuLgYcXFxiIqKgkQigZmZGYqKil7aB5qeKUmlUsHCwqJUCBHhULgYAI7j4OTkhEuXLkGtViM0NBRSqVT7eHx8PBITEwEAdevWhbu7u/axevXqQSqVQiaTYe7cuZg7dy6cnJwQEhKC4ODgUl+Ojo4V/aMRUopKpUJ8fDyio6NLfcXFxWmP0C0tLSGVStGsWTNYWVlpX6tQKHDu3DltiLRs2VJ74STHcfDy8kJmZiZEIhHs7e0F+fnI/6NwMRCNGzfG1q1bkZGRgd9++63UY3PnzsXSpUsBANOnT8fQoUO1j3Ech+3bt8PW1hZ///03srKytA0bERGBH374Qdu0rq6uCA4Ofi547OzsKuznJFWDWq1GQkLCcyESGxurXb7IyckJwcHBaN++PT799FMEBwcjPT0dw4cPx6VLl7BixYpS75mYmIgmTZogLy8Prq6u2LNnT6m/XY7jMGvWLLi5udEpMQNA4WIAMjMzceHCBSgUCmzduhUtWrQotZSFSCQq9d9isVj7/3K5HNu2bUPr1q3Rpk0biMVi9O/fX/u45nRDyQY/duwYNmzYoF2Xyd3d/bnQCQoKKjVYSsiL8DyP+/fvPxcid+7c0Z6asre3R3BwMFq2bIlPPvlE+zf2otO32dnZ8PHxwdatWzF48OBSf4Ml/+45jivVC4wxpKWl4ffff0fdunXpb9cAULgIKCsrCytXrsSGDRsgEonQrFkz7N27F++//z66dev22mmYPM/j119/xfXr13Hw4MFSzadhamqKkJCQ567gLywsRGxsbKkNwqFDhxAWFgb2v0UbPD09Xxg6JU9VkKqB53k8ePDghSGimQJva2uL4OBgNGnSBCNGjND+zbi5ub3xlGJHR0d8+umnmDZtGtatW4evvvrqjdYMKyoqwsKFC5GSkoLExES0aNECCxYsQJcuXd74exPdonARQHZ2NlavXo3169cDAKZOnYrPP/8cxcXF6NWrF8aPH49ffvkFoaGhEIlEEIlEkEgk4DgOHMeBMQa1Wo3du3dj3rx5GD9+PFq1alWmGszNzVG/fn3Ur1+/1L/L5XLcvXu31AZk//79WLVqlfY5tWrVei506tSpAwsLi/L/coigGGNISUl5LkRiYmIgk8kAANbW1ggKCkL9+vUxZMgQ7d/Am16f9TojRozA6dOnsWLFClhaWmL8+PEwNzcHAEgkEkgkklJHLPn5+ViyZAn27NmD77//Hl5eXpg7dy7ee+89bch07NiRQqaC0dpiFSg3Nxdr1qzBt99+C57nMWnSJEybNg1OTk7a58TExGDYsGFISkrC+PHjMXLkSPA8j7S0NACAt7c38vLy8N1332HXrl348MMP8c0338DS0lKvtctkMty5c+e5jc6DBw8APD1N4ePj81zoBAQEaDcMxHAwxpCenv7c5xkdHY0nT54AeDqwXqdOnefG6Ly8vPS+oX78+DEmTpyIw4cPo0uXLpg6dSrq1KmD2NhY8DwPU1NT+Pn54dKlS1i1ahVu3LiBhQsXYvz48RCLxWCM4fjx45g7dy4uXbqE1q1bY+HChQgNDdVr3eT/UbhUAKlUirCwMKxduxZKpRKffvopvvjii5eu4JqamopFixZhz549kEgkCAoKgqenJ9RqNZKSkhAbGwtHR0fMmDEDw4YNg5mZWQX/RP/vyZMniImJeW4DlZqaCuDpGJGfn592w6TZUNWuXRumpqaC1V1VMMaQmZn5whDJzc0F8PQoNjAw8LkQqVWrVqnxvopWUFCAzZs3Y926dXj06BF8fHzg7+8PGxsb5ObmIjY2FmlpaWjcuDHmzZuHdu3aPVcvYwxHjx7F3LlzcfXqVbRr1w4LFixAu3btBPqpqg4KFz168uQJvv32W6xZswaFhYWYMGECvvzyyzdatVWtVuPOnTv4+++/cenSJWRmZsLExATe3t4IDQ1F586dDfp6FqlU+sINWkZGBoCnpzf8/f2fCx0/Pz+YmJgIXL1xKjlTUPN1+/ZtZGdnA3g6/hYQEPBciPj4+LxwvM5QZGRk4MSJEzh16hQSExNRWFgIe3t7hISEoHPnzmjWrNlrj9wZYzh8+DDmzZuH69ev491338WCBQvQunXrCvopqh4KFz3Iz8/H+vXrsWrVKsjlcowbNw4zZsxA9erV3+r9NGMsHMcZ9EbgTWRnZ78wdB4/fgwAMDExQUBAwHOh4+vra/Q/u67k5uY+FyDR0dHaVR0kEkmp36Hmy8/Pz+hvqKVWq8EY045FlhVjDH/++SfmzZuHqKgodOrUCQsWLECLFi30UG3VRuGiQzKZDBs3bsTKlSuRn5+PTz75BF999RVq1KghdGkG72WnbnJycgAAZmZmCAwMfC50vL29BT11o095eXmlTjlqQiQ9PR3A06m5JY/+NF/+/v50yvE1eJ7HH3/8gXnz5iE6Ohpdu3bFggUL0LRpU6FLqzQoXHRALpfju+++w4oVK5CXl4dRo0Zh1qxZ8PT0FLo0o8YYw6NHj7Qb1ZJfeXl5AAALCwvUqVPnudDx8vIymtCRyWQvDJGUlBQATydLlBy3KjlZQsjxtsqA53ns27cP8+fPx507d9C9e3csWLAAjRs3Fro0o0fhUg4KhQKbNm3C8uXLkZOTg5EjR2L27NmoWbOm0KVVapoL5l4UOprpslZWVggKCiq1MQ4JCYGHh4dgU1LlcnmpGXea+pOTk7XP0cy4K/kVGBhI07z1TK1WY+/evViwYAFiY2PRq1cvzJ8/Hw0bNhS6NKNF4fIWCgsL8eOPP2LZsmV4/PgxPvroI8yZMwfe3t5Cl1alMcbw8OHDUhtuzTUaJS/0e1HoVK9eXWehU1hYWOpaIU0t9+/f116gWrNmzedCpE6dOnSBqsDUajV27dqFBQsW4N69e3j//fcxf/581KtXT+jSjA6FSxkUFRXhp59+wtKlS5GRkYFhw4bh66+/hq+vr9ClkVfgeR7JycnPhc6dO3e0K+va2dk9t7EPDg6Gq6vrS0OnqKhIu7ROyfdNSEjQLq3j4eHx3HsGBQXBxsamwn5+UnYqlQq//fYbFi5ciMTERPTv3x/z5s17bqUL8nIULm+guLgYW7ZswZIlS5CWloYhQ4bg66+/Ru3atYUujZSDWq3WrotVMhzu3r2rXVzR0dERQUFBcHd3h5WVFVQqFXJycnDv3j3Ex8drb15VvXr1F4YILQpq3JRKJbZv345FixYhOTkZAwYMwLx581CnTh2hSzN4FC6voFQq8euvv2Lx4sV4+PAhBg8ejK+//hqBgYFCl0b0QKVSISEhAVFRUThz5gyuXr2K+Ph4ZGVl4dk2sbCwgKenJ0JCQtCyZUs0bdoUwcHBcHBwEKh6ok/FxcXYunVrqW3B3LlzERAQIHRpBovC5QU0eyuLFy9GUlISBgwYgLlz59JdHiuJNz1ieXZcpnbt2qWuMdG89nVHMMHBwahWrZqQPzLRkWfPYgwdOhRff/01/P39hS7N4FC4lKA5z7po0SIkJCTQeVYj97ZjLSEhIXBxcXnjAf6SYy8lvxeNvVReJcdfHz16pB1/9fHxEbo0g0Hhgv+fIbJw4ULEx8fj/fffx7x5855bMZgYJs0ssWenJt+5cwcFBQUAABsbmxeGiC5niT1LoVCUuq0BzRqrfJ6dOTpixAjMmTMHtWrVEro0wVXpcNHMbV+4cCHu3r1Lc9sN3Muub4mJiUF+fj6A569v0VxUKeT1Lc963fUuHMfB29ubrncxIs9e8/bxxx9j9uzZ8PLyEro0wVTJcNFclbtgwQLExMSge/fumD9/Pt555x2hSyOoOlfmPys/P/+52xrcvn1be6W+SCSCr68vXalvwAoKCvD9999rV+sYPXo0Zs2aVSVvu1ylwkWzntD8+fNx+/ZtdO3aFfPnz0ezZs2ELq3Kevz48QtD5HVritWqVavKLGT57BpjmtB53RpjtWvXphWmBaJZZ/Cbb76BTCbDJ598gpkzZ8Ld3V3o0ipMlQgXWglVeHl5ebh582aZV0P28fEx+pV89eXZ1ZE1oVNydeTatWuXWmI/JCSEps9WIF2vkG5MqkS4HD58GD179qR7OAho8+bN+OSTT+g+LhXgVfd1qVGjhvY0G6k4eXl5WLduHVavXg0/Pz9cuXJF6JL0rkqECyGEkIplnCOfhBBCDJpBnMxmjCE+Pl57O1ZjJRKJEBISYpTXKNBnIDz6DIRHn4HuGMRpMZ7nMWHCBHh6esLa2lrocspEqVSiuLgYlpaWOHPmDObNm2eUy3ML+RmoVCoA0MnA/enTp+kzEIBSqYRCoYCNjQ31QTm+t1qthkQiKfc1WYbQBwZx5AI8nXI6evRouLq6Cl3KG2OM4dChQ5g7dy66deuGkJCQ5xY4NCYV9RkwxqBWq3H58mXs3r0bUVFR4DgOjRo1wuDBg9GgQYO3mmbMGINMJqPPoIIxxhAeHo4vv/wSffv2RXBwMH0Gb4gxhoyMDPz++++IiIhAbm4uvLy80K9fP3Tt2hUWFhZlDhpD6QODCRdjdezYMdy8eRM+Pj5o0KCB0OUYPMYYsrKysGDBAmzdulV750gAiIiIwM8//4zJkyfjiy++gJWVlcFcVU9ejud57N27F7dv34avry+tcPGG1Go1Dh06hJkzZyI2NrZUGPz+++/o1asXwsLCUKNGDaPsAwqXcpDJZIiMjAQAvPvuu9qrx8mLMcaQlJSEUaNG4eTJk5BIJOjQoQO6d+8OtVqNP//8E+fPn8eSJUuQnp6ONWvWGN3poaooIyMD//77L0QiEfr371/qts3kxZRKJTZt2oTZs2cjPz8fNWvWxIABA+Dj44OzZ8/i4MGD2L9/P7KysrBr1y64ubkZXcDQbLFyiI+Px71792BlZVWprp0pLCyEXC7X6WE1YwzJyckYMmQIIiIi4ODggPXr1+Ovv/7ClClTMG3aNBw5cgRz5syBiYkJtmzZguXLl2vHY4hhYozhv//+Q1paGmrUqIF3331X6JIMnlqtxo8//ogZM2agoKAAvXr1wokTJ7BixQqMHTsWv/76K3755Rc4OTnh1KlTmDZtmnYVb2NC4fKWGGM4deoUCgoKULt2bfj5+Qldkk5kZGRgzJgxmDZtGoqKinTynowxPH78GKNHj8aFCxfg5uaGbdu2YcyYMbC0tATHceA4DjY2Npg9eza++uorcByHtWvX4vDhw4KfOyYvp1Kp8Pvvv4PneXTs2NGoxopeJzMzE6mpqTrfyfrjjz8wc+ZMFBYWYsiQIdi6dSt8fHy0fSCRSNCvXz+EhYXBwsIC+/btw/bt242uDyhc3pJKpcLx48cBAO3atTPKaZcvEhkZiX379uHnn3/Gd999p70fSXnI5XJMmzYN//33HxwcHLB582a89957L1xg0sTEBNOnT0fv3r0hl8sxZ84cZGRklLsGoh+JiYmIjIyEiYkJPvjgA6NdNPRZN2/eRM+ePTFixAhIpVKdvCdjDFevXsWUKVOQn5+Pnj17Yt26dahWrdpzp7xEIhEGDBiAUaNGQalUYvny5Xjw4IFO6qgoleMvQQCpqam4evUqJBIJunTpInQ5OtOzZ0+MGzcOarUaixcvxqlTp8q1x6RSqbBmzRrs2rUL5ubmWLlyJbp16/bK88eWlpZYsmQJPDw8EBMTg40bN+ok5IhuMcZw5MgR5OTkwM/PD82bNxe6JJ3Jzc1FfHw8Tpw4gWXLlpX79CxjDJmZmZg0aRJSU1PRqFEjbNiwAXZ2di/tBYlEgi+//BK+vr64f/8+Nm3aZFR9QOHyFhhjOHv2LLKzs+Hh4YHGjRsLXZLOmJqaYu7cuQgNDUVubi4+//xzpKenv1XAaBYM/eabbwAAU6dOxbBhw167d8txHGrXro2pU6eC4zj89NNPiIuLM7rTApWdQqHA/v37AQA9evSAnZ2dsAXpUNu2bTF79myIxWJ89913CA8PL9ffX3FxMebOnYuLFy/C1dUVGzdufO09hjiOg4eHByZNmgSRSIStW7ciKSnprWuoaBQub4HneRw+fBg8z6NNmzZwdHQUuiSdsrOzw5o1a1CjRg3cuHEDCxYsgFKpLNN7MMYQFRWFqVOnQiaToVevXvjqq6/eeHFKjuMwfPhwBAcH49GjR/jhhx8oXAzMrVu3cP36dVhaWqJv375Cl6NTIpEI48aNQ48ePVBQUIBZs2YhIyPjrXeydu3aha1bt8LExASLFi1C06ZN32j2F8dxGDx4MPz9/ZGeno7ffvvNaPqAwuUtZGRk4OzZsxCLxejRo4fRTRF8HY7jULduXSxYsACmpqbYtm0b9u/f/8Z/1JpTABMmTMDDhw9Rr149rF27tszTih0dHTFx4kSIRCLs2bPHqPbaKjvGGPbv3w+5XI6GDRuiXr16la4PLC0tsXTpUnh4eOD27dtYvnw51Gp1md5Ds5M1e/ZsFBUV4cMPP8Tw4cPLNDbl7OyMjz76CACwc+dOo1mahsKljBhjOH36NNLT01G9enW0bt260jUV8DRghg4dioEDB6KwsBCzZs1641NTcrkc06dPx/nz5+Hi4oKNGzfCy8urzL8njuPQt29f7V7brl27jGavrbLLycnBX3/9BY7j0L9//0p5+2WO4xAYGIhZs2ZBLBZjy5YtOHHiRJn+BqVSKT7//HOkpaWhQYMGWLRoEUxNTctcx4ABA+Di4oL4+Hj8999/RtEHFC5lpFarceDAAfA8j3bt2lWqqZfPMjMzw6JFi1CnTh0kJSVh2rRp2nvVv0xxcTGWL1+uHcBfsWIFWrZs+dYB7OTkhA8//BAAsGvXLp3N3CFvTzPmmJCQAEdHx0p59K6hOT3btWtXyGQyzJ49G1lZWW/0WqVSiWXLluHkyZOwt7fH2rVrUb169bf6XdWqVQudO3eGWq3Grl27jOL6LwqXMnr48CFOnz4NsViMfv36VZqply/CcRy8vLywevVq2Nra4ujRo1i8eDGKi4tf+HylUokNGzZg5cqVAIBp06Zh6NCh5fodcRyHDz74AI6OjoiNjS337DVSfiU3cKGhoahVq5bQJemVpaUlFi9eDDc3N1y7dg0rV6587cadMYY9e/Zgw4YNEIlEmD17Ntq2bfvWISwSiTBo0CCYmJjgzJkzRnGKuPJuGfVAs0Df48ePUatWrUp7SqwkjuPQuXNnfPXVVxCJRFi3bh3Wr19faoCfMQaFQoFVq1Zhzpw5UCqVGDlyZJkG8F/F19cXoaGhUKlU2LNnj1FNx6yMkpKS8N9//0EikWDw4MFvtcioMdGMQc6YMQMikQjff/89jh49+tKdHMYYIiMj8cUXX0ChUOCDDz7A+PHjy72T1aJFC/j4+CAnJwfHjx83+J0sCpcyKCoqwt69e8EYQ/fu3SvdLLGXEYvF+OyzzzBy5EgolUrMmTMHc+fORXp6OuRyOW7duoVRo0Zh3rx5KCoqwrBhw7Bq1SpYWlrq7PsPHDgQYrEYJ0+epNv0CogxhoMHDyIrKwt+fn7l2hs3JiKRCGPGjEH37t0hk8kwdepU3Llz57kNvGYAf9SoUcjIyECzZs2watUqnYxJ2dvbo3PnztrV2Ms6g7OiUbi8IcYYbt68icuXL8PCwgIDBgyoEk2lYWFhgZUrV2LEiBFQqVRYsWIFmjZtihYtWqBt27bYtWsXJBIJpk6divXr18PW1lZnvx+O49CmTRt4enoiMzMTJ0+eNPi9tspKJpNhz549YIyhb9++cHBwELqkCmNpaYnVq1cjICAACQkJ+Pjjj5GYmKj9W+R5HpGRkRg4cCDi4uJQu3ZtbN68+a3HWV6kR48eMDExwZUrVwz+in0KlzfEGMPOnTtRUFCAxo0bo2HDhlUqXDiOg62tLdatW4fly5fDw8MD6enpuHXrFgoLC9G0aVP89ttvWLZsGaytrXX+u3F2dkZoaCh4nsdff/1V5imhpPwYYzh//jyioqJga2uLAQMGCF1SheI4Dr6+vti0aRPc3Nxw8eJF9O7dGzt27MDp06cxZ84c9OnTB7GxsfDz88O2bdsQEhKi052sxo0bo2bNmsjJycGZM2cMeieLltx/Q2lpaTh48CA4jsOHH35YKadevg7HcbCyssLUqVMxaNAgREVFQSqVolatWggJCdFLqJT83j169MC2bdtw4cIFZGZmwt3dXS/fi7yYWq3Gtm3bUFxcjM6dO6NOnTpVagcLePp32K5dO/z6668YO3YsoqOjMWLECIjFYiiVSohEIrRp0wYbNmxA3bp1df77sbe3R+vWrXHv3j0cO3YMw4cPN9gxLwqXN6A5z/zw4UN4eXmhZ8+eVa6pShKJRKhRowZq1KhRYd+T4zg0a9YMrq6uePToES5fvoxevXpV6c+hoiUkJODYsWOQSCT46KOPdDJZwxhpJrmEh4djzZo1iIiIgFwuh5eXF4YMGYIPP/zwlWuGlfd7d+nSRbuTlZOTA2dnZ51/H12gcHkDT548wS+//ALGGD744AO4ubkJXVKV5OrqinfeeQd//fUX/vnnH/Tq1UvokqoMzRImWVlZCAkJQYcOHap0sHMch4CAAGzatAlPnjxBcXExbG1tYWZmptffi2Yny8HBAWlpabh16xZCQ0MN8rOgMZfXYIzh2LFjiIqKgr29PYYPH26QH2RVIBaL0bFjRwDA2bNnUVBQIHBFVYfmjoialRsq0yKVb4vjOIhEItjZ2cHFxQXm5uYVsm1wd3dHSEgIlEolTp48qffv97YoXF5DoVDgu+++g0qlQs+ePavkeWZDwXEcWrVqBUtLSyQkJOD+/ftCl1QlaFa3vnfvHtzc3DBw4EDqAQGZmpqibdu2AIAzZ84Y7JRkCpdXYIzh33//xfnz52FtbY1x48YZ7OBZVeHn54datWqhoKAAly5dMujZMpVFfn4+Nm/eDJ7n0a9fP9SsWVPokqo0zdR8ExMTxMTEGOzN9ChcXkEul2Pt2rUoLi5G9+7d8c4779Aem8Csra3RpEkTMMYMfipmZcAYw9GjR3Ht2jXY29tj9OjR1AMGIDg4GC4uLsjOzkZUVJRB9gGFy0swxvD3338jMjISNjY2+OyzzyCR0PwHoWlOjQHA1atXadxFzwoKCrBhwwaoVCr06tULwcHBFC4GwMnJCSEhIVCr1YiMjBS6nBeicHkJqVSKVatWQalU4v3336ejFgOhuZDM0tISycnJePjwodAlVVqaZUYuXrwIW1tbTJw4kU4LGwiJRIKWLVsCAC5evGiQ4y4ULi/AGMOOHTtw9epVODo64vPPP6ejFgPi7e0Nd3d3yGQy3Lx50yBPCVQGUqkUa9asgVKpRL9+/arcqhSGjOM4NG/eHBKJBHfv3sXjx4+FLuk5FC7PYIwhOTkZa9asAc/z+Pjjj3W6hAMpP1tbW4SEhIAxhsuXLwtdTqWk2cG6du0anJycMGXKFDpqMTB16tSBg4MDsrKyEBsbK3Q5z6FweYZKpcI333yDpKQk+Pv7Y/LkyZX6ni3GSCQS4Z133gEAXL9+3SBPCRgzxhiSkpK0O1gjR46ksRYD5OLiAj8/PyiVSly+fNngjuBpq1kCYwwnTpzAtm3bYGJigpkzZ6JGjRrUVAaG4zg0atQIYrEY8fHxyMnJEbqkSkWlUmH58uVISkpC7dq1MXnyZDpqMUCmpqZo3LgxAODy5csGd58jCpf/YYwhMzMTs2bNQkFBAbp27UoXixmwgIAAVKtWDVlZWUZxVz5jobkh3o4dO2BiYoJZs2ZV6BpypGyaNGkCjuMQHR0NmUwmdDmlULj8j0qlwpIlS3Djxg1Ur14dixcvrpIrHxsLV1dXeHp6oqioCLdu3TK4UwLGiDGGlJQUzJw5E3K5HD169Khy9y0yJhzHoV69erC0tERqaqrB3USPwgVPm2rfvn346aefIJFIMGfOHL0sl010x8LCAnXq1AEA3LhxQ9hiKonCwkLMnj0b0dHR8PLywpIlS2Bubi50WeQVvLy8UL16dRQUFCA6OtqgdrKqfLgwxnDjxg18+eWX2vtdjxw5koLFwHEch/r16wMAoqOjaVC/nNRqNTZt2oTdu3fD3NwcS5YsQWBgIPWBgbOxsUFAQAB4nse1a9eELqeUKh0ujDGkpqZi3LhxSElJQf369bF8+XLaWzMCHMchJCQEYrEY9+/fx5MnT4QuyWgxxnDkyBEsWLAAKpUKY8eOpdNhRkIsFqNhw4YAgJs3bxrUHVqrbLgwxpCbm4sJEybg0qVLqF69Or7//nt4eHhQUxkJPz8/WFlZISsrC6mpqUKXY5QYYzh37hzGjx+PvLw8dO3aFfPmzauyNwIzNhzHoUGDBhCJRIiLizOonawqGS6MMTx58gSTJ0/G4cOHYWtri2+//RbNmzenYDEibm5ucHFxgUKhQHx8vNDlGB3GGK5du4YRI0YgNTUVjRo1wnfffae3uygS/QgMDISVlRUePXpkUIP6VS5cNEcsEydOxK5du2BhYYEVK1agb9++1FBGxtraGt7e3mCM4c6dO0KXY1QYY7h48SIGDx6Me/fuITAwEL/++itq1qxJfWBkatSoATc3N8jlcty9e1focrSqVLhoplp+9NFH2LlzJ8zNzbF06VKMHj2aLhIzQmKxGIGBgQCAmJgYg5opY8h4nkd4eDgGDhyI+Ph4BAQEYMeOHbTMkZGytraGv78/GGOIiooSuhwtgwuXvLw8qFQqnW8oeJ7HhQsX0KdPHxw+fBg2NjZYvXo1Jk6cSItSGrGgoCAAwL1791BUVCRwNbrBGENcXBxyc3N12geMMRQWFmLjxo0YOnQoHjx4gPr162PPnj1o1KgRBYuREovFqFu3LgDg1q1bBnOlvkGFi1QqxZAhQzBt2jRkZmbqpLE04ytr165Fr169cO3aNdSoUQO//PILxowZQ8FixDiOQ0BAACQSCVJTUyGVSoUuSSeSkpLw/vvvo1+/frh586ZONhaMMSQmJmLUqFGYNm0apFIpOnXqhP3796NevXoULEZO8xnGxcUZzD2ODCpcwsPDceLECaxfvx6dO3fG/v37oVAo3ipkNHtp4eHh6NGjB2bMmIHs7Gy0bt0af/31F95//306FVYJeHp6wtraGlKpFOnp6UKXoxO3bt1CRkYGIiIi0LVrV6xZswY5OTlv3Qc5OTnYsGEDOnTogJ07d0IikWDy5MnYvXs3fHx8KFiMnGYny9zcHBkZGXj06JHQJQEwsHDp3bs3wsLC4O7ujqioKAwdOhQ9e/bEgQMHtKcIXtVgjDHwPI+srCzs27cPvXv3Rt++fXHmzBnY2tpi1qxZOHjwIN2XohJxdnaGs7MzCgsLkZiYKHQ5OtGjRw/88ccfaNKkCTIzMzFjxgx06NABP//8Mx49egSe51/bB2q1Gg8fPsTGjRsRGhqKKVOm4MGDBwgMDMT27duxcuVK2NvbUx9UEh4eHnB0dER+fr7B9IFBnROysLDA2LFjERoaiuXLl2Pfvn04ceIETp06BV9fX3To0AHt2rVDYGAgnJycYGZmBuDpve4zMzMRHR2NkydP4tSpU0hOToZarYaVlRX69++PGTNmoGHDhnS0UslYWlrCy8sL8fHxlWY6skgkQps2bXDkyBFs2LABmzZtwo0bNzB27FgsXrwYHTp0wLvvvouQkBA4OzvDzMwMjDHI5XKkp6fjxo0bOHHiBM6ePYuMjAwwxuDq6oqPP/4Yn376KapXr06hUsnY2dmhZs2aSElJQUxMjNDlADCwcAGeHuLVrl0bP/74Iz755BP88MMPOHLkCOLi4hAbG4vvv/8eFhYWsLa2hrm5ubapCgoKUFRUBMYYOI6Dq6srunbtitGjR6NJkyYwMTGhhqqEJBIJfH19ceLECcTFxcHX11foknSC4zg4OTlh7ty5GDRoEH788Ufs27cPDx8+xJYtW/DLL7/A3NwcNjY2MDc3B8/zUCgUkMlk2okNYrEYPj4+GDBgAEaMGAE/Pz+6N1ElZWpqioCAAERGRuL27dvw8fERuiTDCxfgaWOZmJigefPmaNq0KZKTk3HixAn8888/iIqKwqNHj5CXl6e9j4dEIoGVlRX8/PzQsGFDdOrUCe3atUONGjUgEokoVCq5gIAAiMVipKWloWbNmkKXo1MikQgBAQFYtWoVpk+fjv/++w9Hjx7F9evXkZ6e/lwf2NrawsPDA02bNkXXrl3RqlUrODk5UQ9UAZrlkB4+fAhPT0+hyzGccNFc3PiiZSfs7OzQr18/vP/++5DJZMjOzkZ2djbkcjk4joONjQ0cHR3h4OAACwsL7d5ZXl5ehf4MhYWFFfr9dO1Vn4Eha9++PXbs2AEfHx/8+eefQpdTLq/6DMzNzdGtWzd07doVBQUFyMrK0vaBSCSCjY0NnJyc4ODgAHNzc22g5ObmVujPQH0gjHbt2mHbtm2oXbs2/vjjD6HLMYxw4TgONWvWxLp164x6TEShUKBatWpCl/FWKsNnEBkZSZ+BAaDPQFjnz583iM+AYwZwWfPrZoEZE47jjPIUBH0GwqPPQHj0Gejw+xtCuBBCCKlcaOoIIYQQnaNwIYQQonNVIlwYY/jzzz/RoEEDcByHTp064dy5c0KXVamp1Wr89ttvCAgIAMdx6NOnD93rXmAKhQJhYWFwc3ODiYkJxowZg6SkJKHLqtTy8/OxZMkS7Qy+qVOnIiMjQ+iyKgarQtRqNdu/fz8LCQlhAFiXLl3YhQsXhC6rUlGr1WzXrl0sMDCQAWA9evRgV65cEbosUkJBQQFbtWoVc3Z2ZiYmJmzs2LEsOTlZ6LIqlfz8fLZs2TLm4ODATE1N2aRJk1hqaqrQZVWoKhUuGmq1mu3du5cFBQUxAKxbt27s8uXLQpdl1DS/0+DgYAaAvffee+zixYtCl0VeQSaTsRUrVjAnJydmamrKJkyYwB4+fCh0WUZNJpOxb775hjk5OTETE5Mq/TutkuGioVKpSu1l9+zZk129elXosoyK5miwbt262qPB8+fPC10WKYOSe9lmZmZVci+7vAoKCtjq1auZi4sLk0gkdDTIqni4aKhUKrZjxw5Wu3ZtBoD16dOH3bhxQ+iyDBrP8+zgwYOsQYMGDADr2LEji4yMFLosUg55eXls8eLFzM7Ojpmbm7MpU6aw9PR0ocsyaHK5nIWFhTE3NzcmFovZ6NGj2f3794UuyyBQuJSgVCrZ1q1bma+vLwPA+vXrx6KiooQuy6DwPM8OHTrEGjduzACw0NBQdvr0aaHLIjoklUrZggULWLVq1ZiFhQWbNm0ae/TokdBlGRSFQsHWr1/PqlevzsRiMRs5ciRLSEgQuiyDQuHyAkqlkm3ZsoV5e3szAGzAgAEsOjpa6LIExfM8O3LkCGvSpAkDwNq0acP+++8/ocsiepSbm8vmzp3LbGxsmKWlJfvyyy/Z48ePhS5LUIWFhey7775jHh4eTCQSseHDh7P4+HihyzJIFC6vUFxczDZv3sxq1qzJOI5jgwcPZnfu3BG6rArF8zw7duwYa968OQPAWrZsyf7991/G87zQpZEKkp2dzWbPns2sra2ZtbU1mzlzJsvKyhK6rApVVFTEfvjhB+bp6ck4jmNDhw5ld+/eFbosg0bh8gaKiorYpk2bmKenJxOJRGzo0KEsNjZW6LL0iud59u+//7KWLVsyAKxZs2bs2LFjFCpV2OPHj9lXX33FrKysmI2NDZszZw7LyckRuiy9Ki4uZj/99JN2B3PQoEEsJiZG6LKMAoVLGRQWFrKNGzcyd3f3Sn1IHBERwdq2bcsAsHfeeYcdOXKEQoVoZWZmsi+++IJZWFgwW1tbNm/ePJabmyt0WTqlVCrZL7/8wnx8fBgA9sEHH7Dbt28LXZZRoXB5CwqFgq1bt047Q6SyDOadPn2ahYaGMgCsYcOG7NChQxQq5KUyMjLY559/zszNzZmdnR1buHAhy8vLE7qsclEqlWzbtm3Mz8+PAWB9+/ZlN2/eFLoso0ThUg5yuZytXbuWubq6MolEwkaPHs2SkpKELqvMIiMjWceOHRkAVr9+fXbw4EEKFfLG0tLS2GeffcbMzMyYvb09W7JkCXvy5InQZZWJSqViv/32m/ZyhN69e7Pr168LXZZRo3DRgRctp/HgwQOhy3qtCxcusC5dujAALCQkhO3fv5+p1WqhyyJGKiUlhX366afM1NSUOTo6suXLl7P8/Hyhy3oltVrNdu/ezerUqUPLFekYhYsOaZbTcHR0ZKampmzixIksJSVF6LKec/nyZdatWzcGgAUFBbG9e/dSqBCdefDgARs3bhwzMTFhTk5ObOXKlaygoEDoskpRq9Xs999/p+WK9IjCRQ+ePHnCli5dyuzt7ZmZmRmbPHkyS0tLE7osdvXqVdazZ08GgAUEBLBdu3YxlUoldFmkkkpKSmJjxoxhEomEubi4sDVr1jC5XC5oTTzPswMHDrB69eoxAKxz587s3LlzgtZUWVG46FFeXh5btGiRdjmNqVOnsoyMjAqv48aNG6xPnz4MAPP392c7duygUCEVJjExkX388cdMLBYzNzc39u233zKFQlGhNfA8z/7880/tckUdOnRgZ8+erdAaqhoKlwqQm5vL5s+fr11OY/r06WVaToPnefb48WN2//599vjx4zcebI+KimL9+vVjAJivry/bunUrUyqVb/tjEFIu9+7dYyNGjGAikYi5u7uzDRs2sMLCwjd+/dv0Ac/z7PDhw9rlitq3b89OnTpVnh+DvCEKlwqUk5PDvv76a+1yGjNmzHjlchq5ubksLCxMu9aZ5svX15eFhYW99NqC6OhoNmDAAAaAeXt7sy1btlCoEIMRFxfHhg0bxkQiEfPw8GDff/89Kyoqeunz36YPeJ5nR48eZU2bNqXligRC4SKA7OxsNmvWLO1yGrNmzWLZ2dmlnhMeHs6srKwYx3GM47hSTaX5NysrKxYeHq59zZ07d9jgwYMZx3GsZs2abPPmzay4uLiifzxC3sjdu3fZkCFDGMdxzMvLi/3444/P/b2WtQ94nmfHjx8vtVzRP//8Q1PrBUDhIqDHjx+zGTNmMEtLS2ZjY8O+/vprlpOTw8LDw5lYLGYikahUMz37JRKJmFgsZj///DP78MMPmUgkYp6enmzTpk2v3BMkxJBER0ezgQMHMo7jmLe3N/v5559ZcXFxmfpAJBKx5cuXs1atWmmXKwoPD6dQERDHGGO6uWEyeVuZmZlYuXIlNm7cCBMTEygUCqhUKpTlo6levTrmzJmDUaNGwczMTI/VEqIft2/fxoIFC7Bv3z54e3sjNTUVSqWyTH3QsGFDLF68GO+99x44jtNjteR1KFwMSEZGBgYNGoRTp06V+bWrVq3CtGnT9FAVIRXr5s2bGD58OKKiosr82rCwMHz22Wd6qIqUFYWLAWGMwd/fHwkJCWV6Hcdx8PHxQXx8PO2tEaNHfVA5ULgYkKysLDg7O5fr9Y6OjjqsiJCKR31QOYiELoD8P5lMVq7X5+fn66gSQoRDfVA5ULgYEGtr63K93sbGRkeVECIc6oPKgcLFgDg6OsLX17fM54s5joOvry8cHBz0VBkhFYf6oHKgcDEgHMdh0qRJb/XayZMn0yAmqRSoDyoHGtA3MFKpFB4eHlAoFOB5/rXPF4lEsLCwQEpKCuzs7PRfICEVgPrA+NGRi4Gxs7PD/v37wXEcRKJXfzwikQgcx+HAgQPUUKRSoT4wfhQuBqhLly74+++/YWFhAY7jnjvM1/ybhYUFjhw5gs6dOwtUKSH6Q31g3ChcDFSXLl2QkpKCsLAw+Pj4lHrMx8cHYWFhSE1NpYYilRr1gfGiMRcjwBhDTk4O8vPzYWNjAwcHBxq0JFUO9YFxoXAhhBCic3RajBBCiM5RuBBCCNE5ChdCCCE6R+FCCCFE5yhcCCGE6ByFCyGEEJ2jcCGEEKJzFC6EEEJ0jsKFEEKIzlG4EEII0TkKF0IIITpH4UIIIUTnKFwIIYToHIULIYQQnfs/6sawbDiIg94AAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 500x400 with 10 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"model.plot(beta=10)"
]
},
{
"cell_type": "markdown",
"id": "64d2573b",
"metadata": {},
"source": [
"Fix the first layer activation to be linear function, and the second layer to be sine functions (caveat: this is quite sensitive to hypreparams)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "e2e78752",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"r2 is 0.8357976675033569\n",
"r2 is not very high, please double check if you are choosing the correct symbolic function.\n",
"saving model version 0.1\n",
"r2 is 0.8300805687904358\n",
"r2 is not very high, please double check if you are choosing the correct symbolic function.\n",
"saving model version 0.2\n",
"r2 is 0.8376883268356323\n",
"r2 is not very high, please double check if you are choosing the correct symbolic function.\n",
"saving model version 0.3\n",
"r2 is 0.8372848629951477\n",
"r2 is not very high, please double check if you are choosing the correct symbolic function.\n",
"saving model version 0.4\n"
]
},
{
"data": {
"text/plain": [
"tensor(0.8373)"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"model.fix_symbolic(0,0,0,'x')\n",
"model.fix_symbolic(0,0,1,'x')\n",
"model.fix_symbolic(0,1,0,'x')\n",
"model.fix_symbolic(0,1,1,'x')"
]
},
{
"cell_type": "markdown",
"id": "3fae3f32",
"metadata": {},
"source": [
"After setting all to be symbolic, we further train the model (affine parameters are still trainable). The model can now reach machine precision!"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "308b72af",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"pde loss: 1.71e+01 | bc loss: 1.14e-02 | l2: 1.37e-01 : 50%|███▌ | 10/20 [00:11<00:11, 1.20s/it]\n"
]
},
{
"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/var/folders/6j/b6y80djd4nb5hl73rv3sv8y80000gn/T/ipykernel_75424/3364925475.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mtrain\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[0;32m/var/folders/6j/b6y80djd4nb5hl73rv3sv8y80000gn/T/ipykernel_75424/2545871995.py\u001b[0m in \u001b[0;36mtrain\u001b[0;34m()\u001b[0m\n\u001b[1;32m 76\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mupdate_grid_from_samples\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx_i\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 77\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 78\u001b[0;31m \u001b[0moptimizer\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstep\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mclosure\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 79\u001b[0m \u001b[0msol\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msol_fun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx_i\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 80\u001b[0m \u001b[0mloss\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0malpha\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mpde_loss\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mbc_loss\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/opt/anaconda3/lib/python3.9/site-packages/torch/optim/optimizer.py\u001b[0m in \u001b[0;36mwrapper\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m 383\u001b[0m )\n\u001b[1;32m 384\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 385\u001b[0;31m \u001b[0mout\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\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 386\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_optimizer_step_code\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 387\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/opt/anaconda3/lib/python3.9/site-packages/torch/utils/_contextlib.py\u001b[0m in \u001b[0;36mdecorate_context\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m 113\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mdecorate_context\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\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 114\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mctx_factory\u001b[0m\u001b[0;34m(\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--> 115\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\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 116\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 117\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mdecorate_context\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/Desktop/2022/research/code/pykan/kan/LBFGS.py\u001b[0m in \u001b[0;36mstep\u001b[0;34m(self, closure)\u001b[0m\n\u001b[1;32m 441\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mobj_func\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0md\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 442\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_directional_evaluate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mclosure\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0md\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 443\u001b[0;31m loss, flat_grad, t, ls_func_evals = _strong_wolfe(\n\u001b[0m\u001b[1;32m 444\u001b[0m obj_func, x_init, t, d, loss, flat_grad, gtd)\n\u001b[1;32m 445\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_add_grad\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0md\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/Desktop/2022/research/code/pykan/kan/LBFGS.py\u001b[0m in \u001b[0;36m_strong_wolfe\u001b[0;34m(obj_func, x, t, d, f, g, gtd, c1, c2, tolerance_change, max_ls)\u001b[0m\n\u001b[1;32m 48\u001b[0m \u001b[0mg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mclone\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmemory_format\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtorch\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcontiguous_format\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 49\u001b[0m \u001b[0;31m# evaluate objective and gradient using initial step\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 50\u001b[0;31m \u001b[0mf_new\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mg_new\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mobj_func\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0md\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 51\u001b[0m \u001b[0mls_func_evals\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[1;32m 52\u001b[0m \u001b[0mgtd_new\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mg_new\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0md\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/Desktop/2022/research/code/pykan/kan/LBFGS.py\u001b[0m in \u001b[0;36mobj_func\u001b[0;34m(x, t, d)\u001b[0m\n\u001b[1;32m 440\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 441\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mobj_func\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0md\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--> 442\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_directional_evaluate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mclosure\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0md\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 443\u001b[0m loss, flat_grad, t, ls_func_evals = _strong_wolfe(\n\u001b[1;32m 444\u001b[0m obj_func, x_init, t, d, loss, flat_grad, gtd)\n",
"\u001b[0;32m~/Desktop/2022/research/code/pykan/kan/LBFGS.py\u001b[0m in \u001b[0;36m_directional_evaluate\u001b[0;34m(self, closure, x, t, d)\u001b[0m\n\u001b[1;32m 289\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_directional_evaluate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mclosure\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0md\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 290\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_add_grad\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0md\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 291\u001b[0;31m \u001b[0mloss\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfloat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mclosure\u001b[0m\u001b[0;34m(\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 292\u001b[0m \u001b[0mflat_grad\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_gather_flat_grad\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 293\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_set_param\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\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.9/site-packages/torch/utils/_contextlib.py\u001b[0m in \u001b[0;36mdecorate_context\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m 113\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mdecorate_context\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\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 114\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mctx_factory\u001b[0m\u001b[0;34m(\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--> 115\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\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 116\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 117\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mdecorate_context\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/var/folders/6j/b6y80djd4nb5hl73rv3sv8y80000gn/T/ipykernel_75424/2545871995.py\u001b[0m in \u001b[0;36mclosure\u001b[0;34m()\u001b[0m\n\u001b[1;32m 70\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 71\u001b[0m \u001b[0mloss\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0malpha\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mpde_loss\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mbc_loss\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 72\u001b[0;31m \u001b[0mloss\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbackward\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 73\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mloss\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 74\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/opt/anaconda3/lib/python3.9/site-packages/torch/_tensor.py\u001b[0m in \u001b[0;36mbackward\u001b[0;34m(self, gradient, retain_graph, create_graph, inputs)\u001b[0m\n\u001b[1;32m 520\u001b[0m \u001b[0minputs\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0minputs\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 521\u001b[0m )\n\u001b[0;32m--> 522\u001b[0;31m torch.autograd.backward(\n\u001b[0m\u001b[1;32m 523\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgradient\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mretain_graph\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcreate_graph\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minputs\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0minputs\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 524\u001b[0m )\n",
"\u001b[0;32m~/opt/anaconda3/lib/python3.9/site-packages/torch/autograd/__init__.py\u001b[0m in \u001b[0;36mbackward\u001b[0;34m(tensors, grad_tensors, retain_graph, create_graph, grad_variables, inputs)\u001b[0m\n\u001b[1;32m 264\u001b[0m \u001b[0;31m# some Python versions print out the first line of a multi-line function\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 265\u001b[0m \u001b[0;31m# calls in the traceback and some print out the last line\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 266\u001b[0;31m Variable._execution_engine.run_backward( # Calls into the C++ engine to run the backward pass\n\u001b[0m\u001b[1;32m 267\u001b[0m \u001b[0mtensors\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 268\u001b[0m \u001b[0mgrad_tensors_\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mKeyboardInterrupt\u001b[0m: "
]
}
],
"source": [
"train()"
]
},
{
"cell_type": "markdown",
"id": "35985ae9",
"metadata": {},
"source": [
"Print out the symbolic formula"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "f0ec310e",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle - 0.5 \\sin{\\left(3.141592 x_{1} + 3.141593 x_{2} - 4.712389 \\right)} + 0.5 \\sin{\\left(3.141593 x_{1} - 3.141592 x_{2} + 1.570797 \\right)}$"
],
"text/plain": [
"-0.5*sin(3.141592*x_1 + 3.141593*x_2 - 4.712389) + 0.5*sin(3.141593*x_1 - 3.141592*x_2 + 1.570797)"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"formula = model.symbolic_formula()[0][0]\n",
"ex_round(formula,6)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5c3e90ca",
"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.9.7"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Reference in New Issue View Git Blame Copy Permalink