Examples-linear-regression-and-non-linear-pendulum.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Examples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "from sklearn.datasets import fetch_california_housing\n",
    "from dynopt.models.models import Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 1. Linear Regression with Boston Housing Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ".. _california_housing_dataset:\n",
      "\n",
      "California Housing dataset\n",
      "--------------------------\n",
      "\n",
      "**Data Set Characteristics:**\n",
      "\n",
      "    :Number of Instances: 20640\n",
      "\n",
      "    :Number of Attributes: 8 numeric, predictive attributes and the target\n",
      "\n",
      "    :Attribute Information:\n",
      "        - MedInc        median income in block group\n",
      "        - HouseAge      median house age in block group\n",
      "        - AveRooms      average number of rooms per household\n",
      "        - AveBedrms     average number of bedrooms per household\n",
      "        - Population    block group population\n",
      "        - AveOccup      average number of household members\n",
      "        - Latitude      block group latitude\n",
      "        - Longitude     block group longitude\n",
      "\n",
      "    :Missing Attribute Values: None\n",
      "\n",
      "This dataset was obtained from the StatLib repository.\n",
      "https://www.dcc.fc.up.pt/~ltorgo/Regression/cal_housing.html\n",
      "\n",
      "The target variable is the median house value for California districts,\n",
      "expressed in hundreds of thousands of dollars ($100,000).\n",
      "\n",
      "This dataset was derived from the 1990 U.S. census, using one row per census\n",
      "block group. A block group is the smallest geographical unit for which the U.S.\n",
      "Census Bureau publishes sample data (a block group typically has a population\n",
      "of 600 to 3,000 people).\n",
      "\n",
      "A household is a group of people residing within a home. Since the average\n",
      "number of rooms and bedrooms in this dataset are provided per household, these\n",
      "columns may take surprisingly large values for block groups with few households\n",
      "and many empty houses, such as vacation resorts.\n",
      "\n",
      "It can be downloaded/loaded using the\n",
      ":func:`sklearn.datasets.fetch_california_housing` function.\n",
      "\n",
      ".. topic:: References\n",
      "\n",
      "    - Pace, R. Kelley and Ronald Barry, Sparse Spatial Autoregressions,\n",
      "      Statistics and Probability Letters, 33 (1997) 291-297\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Load dataset\n",
    "data = fetch_california_housing()\n",
    "print(data.DESCR)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['MedInc',\n",
       " 'HouseAge',\n",
       " 'AveRooms',\n",
       " 'AveBedrms',\n",
       " 'Population',\n",
       " 'AveOccup',\n",
       " 'Latitude',\n",
       " 'Longitude']"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.feature_names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   MedInc  HouseAge  AveRooms  AveBedrms  Population  AveOccup  Latitude  \\\n",
      "0  8.3252      41.0  6.984127   1.023810       322.0  2.555556     37.88   \n",
      "1  8.3014      21.0  6.238137   0.971880      2401.0  2.109842     37.86   \n",
      "2  7.2574      52.0  8.288136   1.073446       496.0  2.802260     37.85   \n",
      "3  5.6431      52.0  5.817352   1.073059       558.0  2.547945     37.85   \n",
      "4  3.8462      52.0  6.281853   1.081081       565.0  2.181467     37.85   \n",
      "\n",
      "   Longitude  \n",
      "0    -122.23  \n",
      "1    -122.22  \n",
      "2    -122.24  \n",
      "3    -122.25  \n",
      "4    -122.25  \n"
     ]
    }
   ],
   "source": [
    "# Define data names\n",
    "feature_names = data.feature_names\n",
    "target_name = 'MedHouseVal'  # Median value of homes in $1000s\n",
    "\n",
    "# Construct Pandas dataframes\n",
    "X = pd.DataFrame(data.data, columns=feature_names)\n",
    "y = pd.DataFrame(data.target, columns=[target_name])\n",
    "\n",
    "print(X.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   MedHouseVal\n",
      "0        4.526\n",
      "1        3.585\n",
      "2        3.521\n",
      "3        3.413\n",
      "4        3.422\n"
     ]
    }
   ],
   "source": [
    "print(y.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model(['MedInc', 'HouseAge', 'AveRooms'], ['MedHouseVal'], estimator=LinearRegression())\n"
     ]
    }
   ],
   "source": [
    "# Select desired inputs and outputs\n",
    "x_names = ['MedInc', 'HouseAge', 'AveRooms']\n",
    "y_names = [target_name]\n",
    "\n",
    "# Initialize model\n",
    "model = Model(x_names, y_names)  # by default uses a linear estimator\n",
    "\n",
    "print(model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "The feature names should match those that were passed during fit.\nFeature names unseen at fit time:\n- x0\n- x1\n- x2\nFeature names seen at fit time, yet now missing:\n- AveRooms\n- HouseAge\n- MedInc\n",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[21], line 5\u001b[0m\n\u001b[1;32m      2\u001b[0m model\u001b[38;5;241m.\u001b[39mfit(X, y)\n\u001b[1;32m      4\u001b[0m \u001b[38;5;66;03m# Fit score (R-squared)\u001b[39;00m\n\u001b[0;32m----> 5\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[43mmodel\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mscore\u001b[49m\u001b[43m(\u001b[49m\u001b[43mX\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43my\u001b[49m\u001b[43m)\u001b[49m)\n",
      "File \u001b[0;32m~/dyn-opt/src/dynopt/models/models.py:208\u001b[0m, in \u001b[0;36mModel.score\u001b[0;34m(self, X, y, sample_weight)\u001b[0m\n\u001b[1;32m    204\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"Returns the coefficient of determination R^2 of the\u001b[39;00m\n\u001b[1;32m    205\u001b[0m \u001b[38;5;124;03mprediction.\u001b[39;00m\n\u001b[1;32m    206\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m    207\u001b[0m y_selected \u001b[38;5;241m=\u001b[39m y[\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39my_names]\n\u001b[0;32m--> 208\u001b[0m y_pred \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mpredict\u001b[49m\u001b[43m(\u001b[49m\u001b[43mX\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    209\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m r2_score(y_selected, y_pred, sample_weight\u001b[38;5;241m=\u001b[39msample_weight,\n\u001b[1;32m    210\u001b[0m                 multioutput\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124muniform_average\u001b[39m\u001b[38;5;124m'\u001b[39m)\n",
      "File \u001b[0;32m~/dyn-opt/src/dynopt/models/models.py:194\u001b[0m, in \u001b[0;36mModel.predict\u001b[0;34m(self, X)\u001b[0m\n\u001b[1;32m    191\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39minput_scaler:\n\u001b[1;32m    192\u001b[0m     \u001b[38;5;66;03m# Apply the scaling to input data\u001b[39;00m\n\u001b[1;32m    193\u001b[0m     X \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39minput_scaler\u001b[38;5;241m.\u001b[39mtransform(X\u001b[38;5;241m.\u001b[39mvalues\u001b[38;5;241m.\u001b[39mastype(\u001b[38;5;28mfloat\u001b[39m))\n\u001b[0;32m--> 194\u001b[0m y \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mestimator\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mpredict\u001b[49m\u001b[43m(\u001b[49m\u001b[43mX\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    195\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39moutput_scaler:\n\u001b[1;32m    196\u001b[0m     \u001b[38;5;66;03m# Reverse the scaling on output predictions\u001b[39;00m\n\u001b[1;32m    197\u001b[0m     y \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39moutput_scaler\u001b[38;5;241m.\u001b[39minverse_transform(y)\n",
      "File \u001b[0;32m~/anaconda3/envs/torch/lib/python3.10/site-packages/sklearn/linear_model/_base.py:386\u001b[0m, in \u001b[0;36mLinearModel.predict\u001b[0;34m(self, X)\u001b[0m\n\u001b[1;32m    372\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mpredict\u001b[39m(\u001b[38;5;28mself\u001b[39m, X):\n\u001b[1;32m    373\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m    374\u001b[0m \u001b[38;5;124;03m    Predict using the linear model.\u001b[39;00m\n\u001b[1;32m    375\u001b[0m \n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    384\u001b[0m \u001b[38;5;124;03m        Returns predicted values.\u001b[39;00m\n\u001b[1;32m    385\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[0;32m--> 386\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_decision_function\u001b[49m\u001b[43m(\u001b[49m\u001b[43mX\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/anaconda3/envs/torch/lib/python3.10/site-packages/sklearn/linear_model/_base.py:369\u001b[0m, in \u001b[0;36mLinearModel._decision_function\u001b[0;34m(self, X)\u001b[0m\n\u001b[1;32m    366\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_decision_function\u001b[39m(\u001b[38;5;28mself\u001b[39m, X):\n\u001b[1;32m    367\u001b[0m     check_is_fitted(\u001b[38;5;28mself\u001b[39m)\n\u001b[0;32m--> 369\u001b[0m     X \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_validate_data\u001b[49m\u001b[43m(\u001b[49m\u001b[43mX\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43maccept_sparse\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m[\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mcsr\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mcsc\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mcoo\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mreset\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mFalse\u001b[39;49;00m\u001b[43m)\u001b[49m\n\u001b[1;32m    370\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m safe_sparse_dot(X, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mcoef_\u001b[38;5;241m.\u001b[39mT, dense_output\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m) \u001b[38;5;241m+\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mintercept_\n",
      "File \u001b[0;32m~/anaconda3/envs/torch/lib/python3.10/site-packages/sklearn/base.py:580\u001b[0m, in \u001b[0;36mBaseEstimator._validate_data\u001b[0;34m(self, X, y, reset, validate_separately, cast_to_ndarray, **check_params)\u001b[0m\n\u001b[1;32m    509\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_validate_data\u001b[39m(\n\u001b[1;32m    510\u001b[0m     \u001b[38;5;28mself\u001b[39m,\n\u001b[1;32m    511\u001b[0m     X\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mno_validation\u001b[39m\u001b[38;5;124m\"\u001b[39m,\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    516\u001b[0m     \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mcheck_params,\n\u001b[1;32m    517\u001b[0m ):\n\u001b[1;32m    518\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"Validate input data and set or check the `n_features_in_` attribute.\u001b[39;00m\n\u001b[1;32m    519\u001b[0m \n\u001b[1;32m    520\u001b[0m \u001b[38;5;124;03m    Parameters\u001b[39;00m\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    578\u001b[0m \u001b[38;5;124;03m        validated.\u001b[39;00m\n\u001b[1;32m    579\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[0;32m--> 580\u001b[0m     \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_check_feature_names\u001b[49m\u001b[43m(\u001b[49m\u001b[43mX\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mreset\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mreset\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    582\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m y \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m \u001b[38;5;129;01mand\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_get_tags()[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mrequires_y\u001b[39m\u001b[38;5;124m\"\u001b[39m]:\n\u001b[1;32m    583\u001b[0m         \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\n\u001b[1;32m    584\u001b[0m             \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mThis \u001b[39m\u001b[38;5;132;01m{\u001b[39;00m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;18m__class__\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;18m__name__\u001b[39m\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m estimator \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    585\u001b[0m             \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mrequires y to be passed, but the target y is None.\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    586\u001b[0m         )\n",
      "File \u001b[0;32m~/anaconda3/envs/torch/lib/python3.10/site-packages/sklearn/base.py:507\u001b[0m, in \u001b[0;36mBaseEstimator._check_feature_names\u001b[0;34m(self, X, reset)\u001b[0m\n\u001b[1;32m    502\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m missing_names \u001b[38;5;129;01mand\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m unexpected_names:\n\u001b[1;32m    503\u001b[0m     message \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m (\n\u001b[1;32m    504\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mFeature names must be in the same order as they were in fit.\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    505\u001b[0m     )\n\u001b[0;32m--> 507\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(message)\n",
      "\u001b[0;31mValueError\u001b[0m: The feature names should match those that were passed during fit.\nFeature names unseen at fit time:\n- x0\n- x1\n- x2\nFeature names seen at fit time, yet now missing:\n- AveRooms\n- HouseAge\n- MedInc\n"
     ]
    }
   ],
   "source": [
    "# Fit model\n",
    "model.fit(X, y)\n",
    "\n",
    "# Fit score (R-squared)\n",
    "print(model.score(X, y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "The feature names should match those that were passed during fit.\nFeature names unseen at fit time:\n- x0\n- x1\n- x2\nFeature names seen at fit time, yet now missing:\n- AveRooms\n- HouseAge\n- MedInc\n",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[22], line 3\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[38;5;66;03m# Prediction with fitted model\u001b[39;00m\n\u001b[0;32m----> 3\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[43mmodel\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mpredict\u001b[49m\u001b[43m(\u001b[49m\u001b[43mX\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mhead\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m)\n",
      "File \u001b[0;32m~/dyn-opt/src/dynopt/models/models.py:194\u001b[0m, in \u001b[0;36mModel.predict\u001b[0;34m(self, X)\u001b[0m\n\u001b[1;32m    191\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39minput_scaler:\n\u001b[1;32m    192\u001b[0m     \u001b[38;5;66;03m# Apply the scaling to input data\u001b[39;00m\n\u001b[1;32m    193\u001b[0m     X \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39minput_scaler\u001b[38;5;241m.\u001b[39mtransform(X\u001b[38;5;241m.\u001b[39mvalues\u001b[38;5;241m.\u001b[39mastype(\u001b[38;5;28mfloat\u001b[39m))\n\u001b[0;32m--> 194\u001b[0m y \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mestimator\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mpredict\u001b[49m\u001b[43m(\u001b[49m\u001b[43mX\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    195\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39moutput_scaler:\n\u001b[1;32m    196\u001b[0m     \u001b[38;5;66;03m# Reverse the scaling on output predictions\u001b[39;00m\n\u001b[1;32m    197\u001b[0m     y \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39moutput_scaler\u001b[38;5;241m.\u001b[39minverse_transform(y)\n",
      "File \u001b[0;32m~/anaconda3/envs/torch/lib/python3.10/site-packages/sklearn/linear_model/_base.py:386\u001b[0m, in \u001b[0;36mLinearModel.predict\u001b[0;34m(self, X)\u001b[0m\n\u001b[1;32m    372\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mpredict\u001b[39m(\u001b[38;5;28mself\u001b[39m, X):\n\u001b[1;32m    373\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m    374\u001b[0m \u001b[38;5;124;03m    Predict using the linear model.\u001b[39;00m\n\u001b[1;32m    375\u001b[0m \n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    384\u001b[0m \u001b[38;5;124;03m        Returns predicted values.\u001b[39;00m\n\u001b[1;32m    385\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[0;32m--> 386\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_decision_function\u001b[49m\u001b[43m(\u001b[49m\u001b[43mX\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/anaconda3/envs/torch/lib/python3.10/site-packages/sklearn/linear_model/_base.py:369\u001b[0m, in \u001b[0;36mLinearModel._decision_function\u001b[0;34m(self, X)\u001b[0m\n\u001b[1;32m    366\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_decision_function\u001b[39m(\u001b[38;5;28mself\u001b[39m, X):\n\u001b[1;32m    367\u001b[0m     check_is_fitted(\u001b[38;5;28mself\u001b[39m)\n\u001b[0;32m--> 369\u001b[0m     X \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_validate_data\u001b[49m\u001b[43m(\u001b[49m\u001b[43mX\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43maccept_sparse\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m[\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mcsr\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mcsc\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mcoo\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mreset\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mFalse\u001b[39;49;00m\u001b[43m)\u001b[49m\n\u001b[1;32m    370\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m safe_sparse_dot(X, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mcoef_\u001b[38;5;241m.\u001b[39mT, dense_output\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m) \u001b[38;5;241m+\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mintercept_\n",
      "File \u001b[0;32m~/anaconda3/envs/torch/lib/python3.10/site-packages/sklearn/base.py:580\u001b[0m, in \u001b[0;36mBaseEstimator._validate_data\u001b[0;34m(self, X, y, reset, validate_separately, cast_to_ndarray, **check_params)\u001b[0m\n\u001b[1;32m    509\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_validate_data\u001b[39m(\n\u001b[1;32m    510\u001b[0m     \u001b[38;5;28mself\u001b[39m,\n\u001b[1;32m    511\u001b[0m     X\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mno_validation\u001b[39m\u001b[38;5;124m\"\u001b[39m,\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    516\u001b[0m     \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mcheck_params,\n\u001b[1;32m    517\u001b[0m ):\n\u001b[1;32m    518\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"Validate input data and set or check the `n_features_in_` attribute.\u001b[39;00m\n\u001b[1;32m    519\u001b[0m \n\u001b[1;32m    520\u001b[0m \u001b[38;5;124;03m    Parameters\u001b[39;00m\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    578\u001b[0m \u001b[38;5;124;03m        validated.\u001b[39;00m\n\u001b[1;32m    579\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[0;32m--> 580\u001b[0m     \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_check_feature_names\u001b[49m\u001b[43m(\u001b[49m\u001b[43mX\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mreset\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mreset\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    582\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m y \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m \u001b[38;5;129;01mand\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_get_tags()[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mrequires_y\u001b[39m\u001b[38;5;124m\"\u001b[39m]:\n\u001b[1;32m    583\u001b[0m         \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\n\u001b[1;32m    584\u001b[0m             \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mThis \u001b[39m\u001b[38;5;132;01m{\u001b[39;00m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;18m__class__\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;18m__name__\u001b[39m\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m estimator \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    585\u001b[0m             \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mrequires y to be passed, but the target y is None.\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    586\u001b[0m         )\n",
      "File \u001b[0;32m~/anaconda3/envs/torch/lib/python3.10/site-packages/sklearn/base.py:507\u001b[0m, in \u001b[0;36mBaseEstimator._check_feature_names\u001b[0;34m(self, X, reset)\u001b[0m\n\u001b[1;32m    502\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m missing_names \u001b[38;5;129;01mand\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m unexpected_names:\n\u001b[1;32m    503\u001b[0m     message \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m (\n\u001b[1;32m    504\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mFeature names must be in the same order as they were in fit.\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    505\u001b[0m     )\n\u001b[0;32m--> 507\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(message)\n",
      "\u001b[0;31mValueError\u001b[0m: The feature names should match those that were passed during fit.\nFeature names unseen at fit time:\n- x0\n- x1\n- x2\nFeature names seen at fit time, yet now missing:\n- AveRooms\n- HouseAge\n- MedInc\n"
     ]
    }
   ],
   "source": [
    "# Prediction with fitted model\n",
    "\n",
    "print(model.predict(X.head()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'MEDV': 29.01211141973685}\n"
     ]
    }
   ],
   "source": [
    "# Single point prediction is also supported\n",
    "x = {'LSTAT': 4.98, 'RM': 6.575, 'TAX': 296}\n",
    "print(model.predict(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple Pendulum Example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def pendulum_dydt(t, y):\n",
    "    dydt = np.empty_like(y)\n",
    "    dydt[0] = y[1]\n",
    "    dydt[1] = -y[1] - 5*np.sin(y[0])\n",
    "    return dydt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import odeint\n",
    "\n",
    "# Integrate dy/dt over time\n",
    "t = np.arange(0, 10, 0.1)\n",
    "y0 = [0.1, 0]  # Initial condition\n",
    "y = odeint(pendulum_dydt, y0, t, tfirst=True)\n",
    "assert y.shape == (100, 2), y.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk4AAADvCAYAAAAEhluXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABZG0lEQVR4nO3dd3hUZdrA4d+ZkklvpEMCoXcIVUAQFYIgKjZAEEVd1EVXEXVXdPcT3V2xrMIqIlZQEWFXQcFFIDaQJjX0DkkgJISQ3iaTmfP9cUggJkACMzmT5Lmv61wzOXPKM2/KPHmroqqqihBCCCGEuCyD3gEIIYQQQtQXkjgJIYQQQtSQJE5CCCGEEDUkiZMQQgghRA1J4iSEEEIIUUOSOAkhhBBC1JAkTkIIIYQQNWTSOwA9OBwOTp06hZ+fH4qi6B2OEEIIIXSkqir5+flERUVhMFy6TqlRJk6nTp0iOjpa7zCEEEII4UZOnDhBs2bNLnlMo0yc/Pz8AK2A/P39nX59m83G6tWriY+Px2w2O/364uKk7PUjZa8fKXv9SNnrx5lln5eXR3R0dEV+cCmNMnEqb57z9/d3WeLk7e2Nv7+//CLVMSl7/UjZ60fKXj9S9vpxRdnXpPuOdA4XQgghhKghSZyEEEIIIWrI5YnTnDlziI2NxdPTk549e/Lrr79e9Ni0tDTGjRtHu3btMBgMTJkypdrjvv76azp27IjFYqFjx44sXbrURdELIYQQQpzn0j5OixcvZsqUKcyZM4cBAwbw/vvvM3z4cPbt20dMTEyV461WK6GhobzwwgvMnDmz2mtu3LiRMWPG8Pe//53bb7+dpUuXMnr0aNatW0ffvn1d+XZqLCWrCFXVOwohhBANnd1ux2az6R2GLmw2GyaTiZKSEux2+yWPNZvNGI1Gp9zXpYnTW2+9xUMPPcQf/vAHAGbNmsWqVat47733mDFjRpXjW7Rowb///W8APvnkk2qvOWvWLIYOHcq0adMAmDZtGmvWrGHWrFl8+eWX1Z5jtVqxWq0VX+fl5QFaoTv7By6nyEb8rPUEehjZpu4jvmMEvZoHYjJKq2hdKP9+NtY/JHqSstePlL1+9Cp7VVXJyMio+DxrjFRVJSIigpSUlBp16vb39ycsLKzaY2vz/XNZ4lRaWsq2bdt47rnnKu2Pj49nw4YNV3zdjRs38tRTT1XaN2zYMGbNmnXRc2bMmMFLL71UZf/q1avx9va+4liqcyhXQcHAWavC57+d5PPfTuJtVOkYpDKkqYNI595OXERCQoLeITRaUvb6kbLXT12XvZ+fH0FBQYSEhODh4SGTOV+CqqqUlpZy5swZDh06RH5+fpVjioqKanw9lyVOmZmZ2O12wsPDK+0PDw8nPT39iq+bnp5e62tOmzaNqVOnVnxdPl9DfHy806cjGAE8UFjCnCU/k+XVlDWHz5JdZGNrpsL+fDPvjOnGwDYhTr2nOM9ms5GQkMDQoUNlaHAdk7LXj5S9fvQoe7vdzrFjxwgNDaVJkyZ1ck93VD7bd01XAfH09MRisdC/f/8qzXa1qblz+TxOv38zqqpedWZc22taLBYsFkuV/Waz2SU/6AE+0K2JyogRXTEYTWxLzuathINsOpbFpAU7+OeozoztU7WPl3AeV31vxeVJ2etHyl4/dVn2drsdRVHw9fW97PIgDZnD4QC0nKAm5eDr60tmZiZAle9Vbb53LivxkJAQjEZjlZqgjIyMKjVGtREREeH0a7qS0aDQJzaYzx7syx1xTbE7VJ5bsps3Vh1AlR7kQgghrpA0z9WOs8rLZYmTh4cHPXv2rNLum5CQQP/+/a/4uv369atyzdWrV1/VNeuCh8nAm6O78cSNbQB49+ejTFmciLXs0iMBhBBCCOE+XFrHN3XqVD766CM++eQT9u/fz1NPPUVKSgqPPvoooPU9uu+++yqdk5iYSGJiIgUFBZw5c4bExET27dtX8fqTTz7J6tWree211zhw4ACvvfYaP/zww0XnfHIniqIwdWhbXr+rKyaDwreJp3h5+b7LnyiEEEI0cNOnT0dRFBRFueSAr+qUnxcYGOiS2C7k0sRpzJgxzJo1i5dffpnu3buzdu1aVqxYQfPmzQFtwsuUlJRK58TFxREXF8e2bdtYuHAhcXFxjBgxouL1/v37s2jRIubNm0fXrl2ZP38+ixcvdps5nGpidK9o5t7bE4Avfkvhu12ndI5ICCGE0F+nTp1IS0vj4YcfrthntVr505/+REhICD4+Ptx6662cPHmy0nlpaWm1TraulMs7h0+ePJnJkydX+9r8+fOr7KtJv5+77rqLu+6662pD09WQjuE8dn0r3v35KM99vZvOUQG0CPHROywhhBBCNyaTiYiIiEr7pkyZwvLly1m0aBFNmjTh6aefZuTIkWzZsqXimIiICAICAuokxsbbHd8NPDWkLX1aBFNgLeOxhdspsUl/JyGEEA3TZ599RpMmTSpNSA1w5513Vum2Uy43N5ePP/6YN998kyFDhhAXF8eCBQvYvXs3P/zwQ12EXYUkTjoyGQ38+57uBPt4sPdUHq+s2K93SEIIIeoZVVUpKi3TZavN6PC7774bu93OsmXLKvZlZmby3Xff8cADD1R7zrZt27DZbMTHx1fsi4qKonPnzmzcuPHKC+0quLypTlxaZIAXb43uxsR5W/hsYzJ9Y5twc9dIvcMSQghRTxTb7HT8v1W63Hvfy8Pw9qhZKuHl5cW4ceOYN28ed999NwBffPEFzZo1Y/DgwaxZs6bKOenp6Xh4eBAUFFRp/9VOpn01JHFyA4PbhfHHwa1475ejPPf1Lno0DyQywEvvsOqnMweJyt6MkpgN9mKwFkBZMURfA62HQCOeLE4IIfQ2adIkevfuTWpqKk2bNmXevHlMnDix1nMsOWMy7SsliZObeHpoWzYePUviiRxmrDjA2/fE6R1S/ZJ/Gn58CXPiF/QGSKrmmCatoe+j0O0esPjWbXxCCOEiXmYj+14eptu9ayMuLo5u3brx2WefMWzYMHbv3s3y5csvenxERASlpaVkZ2dXqnXKyMigX79+Vxz31ZDEyU2YjAb+Maozt8xex7KdpxjfN4a+LRvvGkQ1VlYKm9+HX16DUm3hxiyf1gRGtMDg6a8lSA4H7F8GZ4/Aimfgx79Dz/tg4DPgFahv/EIIcZUURalxc5k7+MMf/sDMmTNJTU1lyJAhREdHX/TYnj17YjabSUhIYPTo0YA29cCePXt49dVX6yrkSqTdwo10bhrAPefWsHtx2V7K7A6dI3JzyRth7gBY/VctaYqKo2ziSn5t+3/Yxy6Cu+fBre/AqHdh6j4Y/gYEtwRrLmx4Bz6Oh+wkvd+FEEI0KuPHjyc1NZUPP/yQBx988JLHBgQE8NBDD/H000/z448/smPHDu699166dOnCkCFD6ijiyiRxcjPPxrcjwMvMgfR8Fm5OufwJjVXqNvj8dsg8BN4hcOts+MNPqE17VX+8xQ/6PgyPb4N7FoFfFGQehA9vhBOb6zZ2IYRoxPz9/bnzzjvx9fVl1KhRlz1+5syZjBo1itGjRzNgwAC8vb1Zvnw5RmPtmgmdRRInNxPk48Ez8W0BeHP1IbIKS3WOyA1lJ8HCMVqn71Y3wJ+2QY8JNev4bTBAu+Ew6UeI6ApFmTB/JOxZ4vKwhRBCaNLS0hg/fjwWi+Wyx3p6evLOO+9w9uxZioqKWL58+SWb91xNEic3NK5vczpE+pNbbOONVQf1Dse9FGXBgrug8AxEdIHRn11ZPyX/KHjge2g7HOxW+OoB+PVNqMWcJEIIIWonKyuLRYsW8dNPP/HYY49VeX337t34+voyZ86cWl3X19e3Yh1cV6s/vckaEaNBYfotHRnzwSYWbUlhXJ8YujSrm6nk3VqZFRbfC2cPg39TGPdfrQnuSll8YewXWh+pTXPgx5fB5AX9ql8iSAghxNXp0aMH2dnZvPbaa7Rr167Sa0888QT33nsvAKGhobW6bmJiIkCdNN9J4uSm+rZswq3doli28xTTl+/lq0f76TZnhVtwOOCbyZC8Hiz+MP6/4O+EiUINRrhpBviGwQ/TYfUL2rQFbeMve6oQQojaSUpKuuhrwcHBBAcHX9F1W7dufYUR1Z401bmx50d0wNNsYFtyNj8fzNA7HH2tnwV7vgKDCcZ8DuGdnHv9AVOgx32gOuCrByFDlr8RQghRlSRObiwiwJP7+7UAtI7iDkcj7X+TcwLWvK49v/lNaDnY+fdQFBjxJjS/VpvaYOEYKMx0/n2EEELUa5I4ublHrmuFj4eRvafyWLVXn3V5dJfwN20EXUx/6HG/6+5j8tBqs4JiISdZ609VZr38eUIIIRoNSZzcXLCPBw9dGwvAWwmHsDe2Wqfjv8LepaAYYPhrWs2QK3kHw7jFWj+qlI3w3VTX3k8IIUS9IolTPfDQwJb4e5o4nFHA8p2n9A6n7tjL4Pu/aM97PgCRXevmvqHttFnHFQMkLoDdX9XNfYUQQrg9SZzqgQAvM49c1wqAWT8cajxLsWybBxl7wTMQbvhr3d679RAY9Kz2/H9TIa8RJaxCCCEuShKnemJi/xYE+3iQdLaIJdtT9Q7H9Yqy4Kd/aM9v+KvWhFbXBj0LUXFQkgvfPiaTYwohhAtNnz4dRVFQFIVZs2bV6tzy8wIDA10S24VcnjjNmTOH2NhYPD096dmzJ7/++uslj1+zZg09e/bE09OTli1bMnfu3Eqvz58/v6KALtxKSkpc+TZ052Mx8cdztU7//vEw1jK7zhG52E//gJIcCO+sNdPpwWiG2z8Akycc/Qm2fKRPHEII0Uh06tSJtLQ0Hn744Yp9H3zwAYMHD8bf3x9FUcjJyalyXlpaWq2TrSvl0sRp8eLFTJkyhRdeeIEdO3YwcOBAhg8fTkpK9YvXHj9+nBEjRjBw4EB27NjB888/zxNPPMHXX39d6Th/f3/S0tIqbZ6enq58K27h3muaE+ZnITWnmP9sOaF3OK5zeq/WTAdah3CjjvO0hraFIS9pz1f/DTIP6xeLEEI0cCaTiYiICLy9vSv2FRUVcdNNN/H8889f9LyIiAgCAupmhQ2XfiK99dZbPPTQQ/zhD38AYNasWaxatYr33nuPGTNmVDl+7ty5xMTEVGSNHTp0YOvWrfzrX//izjvvrDhOURQiIiJqHIfVasVqPT+sPC8vDwCbzYbNZruSt3ZJ5dd09rVNCvzxulhe+u4A7/58hDu6R+JhanitrcZfZ2JQHTja34K9aV+oRTm6pOx7PIDxwP8wJK3FseRh7Pev0CbiFJW46udeXJ6UvX70KHubzYaqqjgcDhyO+tPn9bPPPuPpp5/m5MmTlRb3veuuu/Dx8aFly5YAVd7TE088AcAvv/xS8brD4UA9133i92VxsTIpP8dms1VZmqU23z+X/fUvLS1l27ZtPPfcc5X2x8fHs2HDhmrP2bhxI/HxlZe6GDZsGB9//DE2mw2z2QxAQUEBzZs3x2630717d/7+978TFxd30VhmzJjBSy+9VGX/6tWrK2W1zpaQkOD0a/o5wN9sJD3Pyt8/X0W/8IbV78azNIuhe5cA8Kvak5wVK67oOs4ue0+f27nBuBXzqe0cnD+ZQxGjnHr9hsQVP/eiZqTs9VOXZV9eK1NQUEBpaanW/7KsuM7uXzkYrxpPEzNs2DCefPJJFi9ezKhRowA4e/Ys//vf//jqq69Yv349dru9onLj94qKigDIz8/HYDhfaZCfnw9ASUkJqqpe9PzS0lKKi4tZu3YtZWVl1V67JlyWOGVmZmK32wkPD6+0Pzw8nPT06idyTE9Pr/b4srIyMjMziYyMpH379syfP58uXbqQl5fHv//9bwYMGMDOnTtp06ZNtdedNm0aU6een48nLy+P6Oho4uPj8ff3v8p3WpXNZiMhIYGhQ4dWJHvOdDYoiRkrD7Ehx5cX7xuAydhwap0MP72MATuOmH70v+vxWp/vyrJX9njCt3+kfcZ3tL7tWWhS/c9bY+Xqn3txcVL2+tGj7EtKSjhx4gS+vr5aN5XSQgyvdqiTe/+e47mT4OFTo2P9/f0ZN24cixcv5r777gO0fsvNmjVjxIgRbN26FaPReNHP5fKKDj8/P/z9/VFVlfz8fPz8/FAUBU9PTxRFuej5JSUleHl5MWjQoCrdey6WbFXH5e0Nv1+YVlXVSy5WW93xF+6/5ppruOaaaypeHzBgAD169OCdd97h7bffrvaaFoulUrVgObPZ7NIfdFddf0L/WOauPU5KVjGr9mcyKq6p0++hC2sB7PgUAEP/JzBcRdm5pOy73wN7l6AcScC88s9w/3LXT8hZD7n690pcnJS9fuqy7O12O4qiYDAYtJoXg37/PNf2/g8//DC9e/cmLS2Npk2bMn/+fCZOnIjRaKz4nDdc5Hrl+8vfd3mTXKWyuMz5iqJU+72qzffOZYlTSEgIRqOxSu1SRkZGlVqlchEREdUebzKZaNKkSbXnGAwGevfuzeHDjafTrreHiYeujeVfqw/x7s9HuLVbFAZDA/gAT/xCG/of3BLa3qR3NFUpCtz8L3j3Gkj6FXYu0pIpIYTQk9kbntdprjlz7bq7xMXF0a1bNz777DOGDRvG7t27Wb58uYuCcw2XpakeHh707NmzSrtvQkIC/fv3r/acfv36VTl+9erV9OrV66LZoKqqJCYmEhkZ6ZzA64kJ/VrgZ9FmE1+9rwGsYeeww6Y52vNrJuv6H9QlBbWAwedmM1/9gjbflBBC6ElRtOYyPbYrqHX/wx/+wLx58/jkk08YMmQI0dHRLigU13Hpp9PUqVP56KOP+OSTT9i/fz9PPfUUKSkpPProo4DW96i8nRPg0UcfJTk5malTp7J//34++eQTPv74Y5555pmKY1566SVWrVrFsWPHSExM5KGHHiIxMbHimo1FgJeZ+/u3AGD2z0cqmjTrrYMrIDsJvIKg+zi9o7m0fo9DWEcoOqstQCyEEKLGxo8fT2pqKh9++CEPPvjgZY9PT08nMTGRI0eOALB7924SExPJytLnH1eXJk5jxoxh1qxZvPzyy3Tv3p21a9eyYsUKmjdvDmgTVl04p1NsbCwrVqzgl19+qRgt9/bbb1eaiiAnJ4eHH36YDh06EB8fT2pqKmvXrqVPnz6ufCtu6cFrY/EyG9mTmscvh87oHc7V2TBbe+z1YI07GurGaIaRs7TnOxZA0npdwxFCiPrE39+fO++8E19f34rRdZcyd+5c4uLimDRpEgCDBg0iLi6OZcuWuTjS6rm8c/jkyZOZPHlyta/Nnz+/yr7rrruO7du3X/R6M2fOZObMmc4Kr14L9vFgfN8YPlp3nNk/HWFw29BLdrx3Wye3wolNYDBDn4cvf7w7iOmrzWi+bR58NwUeXQemqgMQhBBCVJWWlsb48eOrHbj1e9OnT2f69OlV9jscjlqNhnMWN+1IImpq0qCWeJgMbEvOZtOxetrfZuO52qYud4NfzSc21d2QF8EnFDIPwYbqR3QKIYQ4Lysri0WLFvHTTz/x2GOPVXl99+7d+Pr6MmfOnFpd19fXt8667Mj0x/VcuL8no3s1Y8GmFOb8coR+raoffei28k7Bvm+15/2qr5l0W15BMOwVWDIJ1r4JXcdAYIzeUQkhhNvq0aMH2dnZvPbaa7Rr167Sa0888QT33nsvAKGhobW6bmJiIkCVGcFdQRKnBuCRQa34cvMJfj2cyc4TOXSLDtQ7pJrbtRhUB8T0g4guekdTe13uhm2fQvI6WPU8jFmgd0RCCOG2kpKSLvpacHAwwcHBV3Td1q1bX2FEtSdNdQ1AdLA3t3WPAuDdn4/oHE0tqKo2FxJAt3o6H5KiwIg3QDHC/uVw5Ee9IxJCCOFCkjg1EJMHt0JRYPW+0xw6na93ODVzagecOQAmT+g0Su9orlx4R+j7iPb8+z9DmfXSxwshhBPU+2lo6pizyksSpwaidZgfN3XSOlbPqS+1Tju/1B7b3wyeAfrGcrUGPwc+YXD2CGx8V+9ohBANWPmE0LVZmFacL6+rXRpH+jg1II9d35rv96SzbOcppg5tR0yT2k2FX6fKSmH3V9rzbm4+4WVNeAZA/N9h6SOw9g2t71Ng/ZoNVwhRPxiNRgIDA8nIyAC0xW/r5VQ0V8nhcFBaWkpJSclF16cDraapqKiIjIwMAgMDr7oDuSRODUjnpgFc1zaUNYfOMHftUV653Y07Wx9eDcVZ4BsBLQfrHY1zdB0D2+ZDykZtOZbRnwFgLbNzIC2ftNxiCqx2CkpsFJbaKbSW4etpItTXQpi/J2F+FsL9PQn28dD3fQgh3F5EhNbCUJ48NUaqqlJcXIyXl1eNEsfAwMCKcrsakjg1MI/f0Jo1h87w1daTPHFDGyICPPUOqXrlzXRd7wZjA/kxPNdRXH1/EMq+b/lswTyW5LZl36k8Su2OGl+maaAXPZoH0TMmkJ7Ng2kf6YfZKK3qQojzFEUhMjKSsLAwbDab3uHowmazsXbtWgYNGnTZ5jez2ey0qQoayCeWKNe7RTB9WgSzOSmLj349xl9HdtQ7pKoKz8KhVdrzhtBMB9jsDn46kMF/t1q5tiyeicaVDDj0Gv8ofZVSzAR5m4kN8cHX04yfxYSPxYi3h4n8kjIy8ks4k2/lTL6Vs4WlpOYUk5pTzPKd2mrnvhYTN3YIY0SXSK5rG4qn2fXzlAgh6gej0Vgncxe5I6PRSFlZGZ6enlfdb6k2JHFqgCZf34rN87L44rcUJl/f2v2afvZ8DQ4bRHbTRqTVY6k5xcxbd5xvElPJLCgF4Dfu4lbTb7QypPFN3HZ8hvyFmOCa9UEosJax60QO25Kz2ZaSzfbkbPJKyvg28RTfJp7Cx8PIjR3CubVbFNe3D8NoaHz9GmrN4YCCdMg9CbknIDcV8lK10Y+KAQzGc48mbeb6gGYQEKP1UfMJvaLV34UQDZckTg3QdW1D6dI0gN2puXz06zH+fFN7vUOqbOdC7bG+zt0E5BbZmPPLEeZtSKK0TGuGC/G1cGePptzdqxnB6SosmUTHwx9A/EOgNK/RdX0tJvq3DqF/6xAAHA6VxJM5rNiVxordaZzKLWHZzlMs23mKpoFejOsbw+he0YT6yTp5FRx2SN8NSesgeb22leRe2bU8fCEqDpr1gqa9oFlv8At3brxCiHpFEqcGSFEUnrixDZM+28qnG5KYNLAlQe5S65RxQJu/yWCCznfpHU2tldjsfL4xmdk/HyG3WOtX0Dc2mEkDW3Jdu9DzfZFCL5hRfOU0uGfhFd3PYFDoERNEj5ggXri5A4kncvhuVxpLtp8kNaeYN1YdZNYPh7ipcyQT+7egZ/MgZ73V+sXh0BKkHQvg4Pdg/V2ipBjBP+pcbVIz8G8KZm9t1vryzW6FvDStVionBfLTobQAkn7VtnJBsdB2GLSJhxbXyuLOQjQykjg1UEM6hNEpyp+9p/L4aN0xnh3mJrVO5Z3CWw8F39qtRaS39Ucy+fNXu0jNKQagbbgv04Z3YHC70KrNcIoCN/8L5l4LB/+n9elqO+yq7q8oCnExQcTFBPHssHas2J3G55uS2ZGSw/Kdp1i+8xR9YoP543Wtqo+pIco5of1M7VgAOcnn93v4QfN+WmLT4lqI6ArGWvaBKCvV5uVK3Qont8DJbZCxD7KPw29ztc3sA62uh/Yjz81H5u/c9yeEcDuSODVQiqLw5I1tePjzbXy6IZk/XOsGtU4Ox/m5m7rXn2a6Epud11ce5JP1xwGI8Pdkanxb7uzR7NJ9jMI6wDV/hA3vwIpnIXYQmL2cEpOn2cgdPZpxR49m7EnNZf6GJL5NTGXz8Sw2H8+ifYQfj17XipFdIzE1xBF5mYdhzWvnfp7OzQbs4Qdd7tSagJv2uvrRmiYPrQ9eeEfocZ+2ryQPjq+Fw6vg0Gqt79SB77TN5Kklx13u1v4xMLvpiFYhxFWRxKkBG9oxnI6R/uxLy+Pjdcd5Zli7y5/kSqe2Q95Jrd9Im6urfakre0/lMmVRIoczCgAY3zeGF27ugLdHDX91rnsOdn+t1YasmwnXP+/0GDs3DeBfd3fj6fi2fLLuOAt/S+FAej5TFicy84dDPHZ9a26Pa9owpjQ4e1SbYLR8cWiAFgMhbgJ0uAU8XDzpq6c/dBipbQ4HpO/Smgb3fA1nD8O+b7XNEgCdb4e4+6BpD+lgLkQD0gD+koqLKe/rBDB/QxI5RaX6BrTvW+2xTbzb/zeuqirvrznKqHfXczijgBBfC59M7MU/b+9S86QJwOILN72iPV83E84cck3AQGSAFy/c3JENz93IM/FtCfbxIPlsEX/+ahfX/+sXvtycUtGRvd4pOAPfPg6ze2tNc6oD2o2AR36Fid9BtzGuT5p+z2CAqO5w/TR4fAs8shb6/0nrP2XN1SZD/egGmNMPNszW3oMQot6TxKmBi+8YTodIfwqsZXyy7rh+gagq7F+mPe94q35x1ECJzc4TixKZ8f0BbHaVoR3DWTVlIDe0v8LRVB1HQeshYC+F5U9qNRUuFOBt5vEb2rDuL9fzwogOhPh6cDK7mGlLdnP9v35hwaZkrGV2l8bgNA6HloDM7gk7PgfVriXek36Ge76EyK56R6hRFG16jfh/wJQ9cP9ybSZ5kyec2a/NJP9We/jPfXDkR5f/DAghXEea6ho4g0HhyRtb8+iC7cxbn8RD17YkwLvuJgqrkL4bspPA5KX1/3BTp/NKePizrew8mYvJoPDirZ24t2/M1XW0VhQYORPevQZSNsC2edD7IecFfRHeHiYmDWrJvdc058vNKcxdc5TUnGL++s0eZv90hD8ObsWY3tHuO6Hm6X3w3VNwYpP2dUQXGPEmxPTVN67LMRi0/myxg2D467B3CWz/XGuqPteUV+oXTWrs3RxpehvZhiYUlZZRZLNTUmrH5lBxqCqo4FBVVBXMJgMeRgMeJgMWkwFPsxE/TxN+niZ8LWa8TJBtBavNXqcTAQrRGLk8cZozZw5vvPEGaWlpdOrUiVmzZjFw4MCLHr9mzRqmTp3K3r17iYqK4s9//jOPPvpopWO+/vpr/va3v3H06FFatWrFP//5T26//XZXv5V6K75jBO0j/DiQns/H648zdWjbug+ivLap9Y1a85Ub2nUyh0mfbeV0npVAbzPvje9Jv1ZNnHPxwBi48f9g5V8g4UVoN1wbHl8HvDyMPHhtLOP6xrB4ywne++Uo6XklvLhsL+/+fISHB7VkXN+Y2jVBulJZKax5Fdb/Gxxl2si1G16APo+4/fI8JTY7xzMLOXamkBPZRaRmF5Oa04WThdPx4yC3lCVwh3Ed/vkniN31FtE7Z/GzozsJ9sH87OhO2VX9STYxffuP+HgYCfb1INjHQqivhXB/bQ3E8rUQIwI8iQrwwt/L1DhGXpaVap34C89AYeb5x+IssBaANV+bdsKaD2UlWs2w3XbusRRQtH9+FMP5zWTR/gk0e2I0eNDrbB7G71aBVyBY/LTNMwC8gsA7WHv0Ctae13Z0p3A7Lv0rtHjxYqZMmcKcOXMYMGAA77//PsOHD2ffvn3ExMRUOf748eOMGDGCSZMmsWDBAtavX8/kyZMJDQ3lzjvvBGDjxo2MGTOGv//979x+++0sXbqU0aNHs27dOvr2dfP/RHWi1Tq14Y9fbOeTdceZ2L9F3c8mvq+8me62ur1vDa3YncbU/yRSYnPQJsyXj+7vRfMmPs69SZ9JsPu/2vD2/z0NYxfWaadhT7OR+/u3YGyfaP679STv/aLVQP3jf/uZ/fMRJvZvwf39Wug7+jLzCHz9EKQlal+3HwnDX9PmXnIjdofK8cxC9p7KZd+pPA5nFHAko4AT2UWo6sXOaso2JvKGYxx3eW5jtPIDnez7GWrczlDjdvJMwewMHsHO0FvI8YrBYFAqfjxsZSqldjtWm4NSu4PiUjsF1jLyS8rIL7GRX1JGdpEVh6poC0hnFXMiq/iS78HLbCQyUEuiogI9iQr0IirQi6bnHiMDPN23NvJCDoc26OTsUcg6pm25J89vBaepGHnpAgagKUDO5pqd4BmozUjvEwo+TcA3HHzCwLd8C9cefcLcvi9oY6Wo6sV/za9W37596dGjB++9917Fvg4dOjBq1ChmzJhR5fi//OUvLFu2jP3791fse/TRR9m5cycbN24EYMyYMeTl5fH9999XHHPTTTcRFBTEl19+WaO48vLyCAgIIDc3F39/J8+7oqrYslJYv/obBtw12W2qzR0OlVtmr2PvqTweGNCCF2/pVHc3zzgAc/qCwQx/Pqr9J+YiNpuNFStWMGLEiBqX/ZebU3h+6W5UFa5vF8rb98Th5+mi79vpffD+QK0m5e5PodMo19ynBkrLHCzZfpK5a46SdLYIAG8PI/f0ieHBa2NpGli7qROupOwrqKo2F9P3fwZbkfbhcsu/dS2fcqqqkpJVxI6UHBJP5LDrZA770/IptlXfT8zf00TrMF9igr1pGuRFsyDvc8mIJyG+Fvw9zRjKp7E4cwh2fAaJX0JR5vmLNO0JXcdC5zvAJ+SyMdpsNv73vxUMvGEoeaUqWYVWzhaUcqbAyuk8Kxl5JWTkW0nPLSE9r4SswpoNFGni41GRREUFehER4ElkgCcR/p5EBngR5m+pu+TKbtPm1TpzQPubcuYAnDmoJUp266XPNXqcS1RCzictXkFg8ddqwC1+2mhfs7d2rNF8/hEqT5TqKNNqscqKUW0lFBflsSdxOy2jI3FYC3CU5KGW5KFYczFZczCX5mApzcVSlodSywQuDx/OEshZAsgkkEw1gDOq9jxL9ecsAWThTxb+FOOJQQGDoiXcBoOCUVEwGBRMBgXjBY9mowGz0YDJqGA2GDCbzu/zONcsbDYq554btcdzzcQeRgMW8/mm44rjTQYsF+yrfC3teuVNzuVxXG2N51X9zfl9WdciL3BZjVNpaSnbtm3jueeeq7Q/Pj6eDRs2VHvOxo0biY+Pr7Rv2LBhfPzxx9hsNsxmMxs3buSpp56qcsysWbMuGovVasVqPf+LlZeXB2iF7uxVpZWUDZg/v5XeHqHYbJOceu2r9czQNjzw6TYWbErm3j7NiAmum1FIhj1LMQKO2OuwG73BhSt5l38/a/p9/eDX47yx+jAAY3s3Y/rIDhgNNT+/1oLbYOj/JMZ1b6KueJay6P7aH3AdKMCdcZGM6hbBqr2nmbv2OPvT8/l43XHmrT/O0A5hTLgmhj4tgmr0B662ZV+hOAfj909j2K+NunQ0vxb7rXO0pkwdVn232uzsSs1jW3I220/kkHgil+yiqnF4mQ20i/CjY6QfbcP9aB3qQ6tQH5r4eFyyvOz2MuzlOVdgLFz/IgyahnJ4NYbEBSjHfkZJ3Qap21BXTUNteQOOjqNQ29x00Qk2bTYbigJeJvD38qBZwKVrDUtsdk7nWUnLLeFUbjFpuVbScotJyy0hNaeEtNwSikrtnC0s5WxhKbtTL75kjb+niTA/C2F+FkL9LDTx8SDYx4MmvtpjkLeZAE8zAd5m/D1NNZsWo+gsyuk92paxD+X0Xsg8iOKo/udBNZghqDlqUEvU4FgIiEH1bwb+UagBzcA7pFLtrsOhUmzTau0KreceS8soKLGTb7VRYLWTX1JWUatXYC2joEQl36pSUKJSYFUosHpQYDXgUL2BEZB06bdkwEEgBQQreYQoeQSTRxMljxAll1ByCFVyz205hJCLRSnDn0L8KSSWVO0iyrmtGsWqB9n4kqP6ke3wJdvuR47qQx4+5KneFY8FeFGoelKEJ9lYKFQ9KcGCFTOlmC5+gyumYsGGhVLtUbHhSSlelOJnLMVbseFjtOGrWPFVSvBRSvFWrHhjxQsr3pTgdcGW0+Yu+t7xBHAVf3OqUZtruCxxyszMxG63Ex5eeSRSeHg46enp1Z6Tnp5e7fFlZWVkZmYSGRl50WMudk2AGTNm8NJLL1XZv3r1ary9nZs8eJTlMxzwKT3Ddyu/w250r+UY2gUYOJhr4NnP1nJ/27oZ2TP4wEICgJ225qSsWFEn90xISLjk66oK350w8EOq9kd8SFMH1xiTWLUyyeWxGRwdGWyJxK8wjbT5D7Cj+cMuv2dNPNICDgQp/JiqcDjPwKp9Gazal0Gkt8qgCAc9Q1QsNahcuFzZXyigKInex9/GpzQTB0b2R93JkaARsC4RSLzCd1I7RWVwLF/hWJ7CsXyFlAKwq5U/PIyKSrQPNPdTifFRaeajEuYFBuUscBYy4WwmnN1f/T1qxgD+92HpdBtNszcRnbWewOIklCOrMRxZjUMxcsavE6cCe5Me0INSk1+VK9Sm7Mt5ArFArBkI0TZVhWK71uE826qQXQo5VoWcUrTt3PMyVSGvpIy8kjKOnCms0f0sBhVPI1iM2vNoYyYdlGTaOJJopSYT60gmRD1b7bnFiidpxqacMjbllKEZqcampBoiyVBCKLMZsJ+GsnSwOcq309gcGZQ6oNQO1nOPpQ7nJgcK2u+GxYj23gxgMap4GMHDoO33MICHwRezwQezIRKzAWwGyDRAjgGOK2AygFEBIw681SL8HDn42nPxLt/K8vC25+JVlounPR+vsjw87XmYVBteSileZBGlZF3Ve7FhxqaYsWPErhi1R4zYMaACqqrVm6koKDhQVBUDDhQcGHBgxI5ZLcOEHRNleCg1HMGrUqMW1e+TdrDid58jV/Jz/3tFRUU1PtblPS1//5+XqqqX/G+suuN/v7+215w2bRpTp06t+DovL4/o6Gji4+Od31QHqEf/D6XoLDd2j8YU3cvp178azbvnMeq9TWw/a+Bv3frTuamLl4jIPo55RwqqYqTznX+hs3ewS29ns9lISEhg6NChF626dThUXvrffn5IPQnAs/FteHhgrEvj+j2lWwTqZyOJyVpH1OAHUDu4R9+vm4GngUOn8/n8txN8m3iKtCIHi48ZWXbCwLCO4YyKi+Ka2OAqs6bXpOwvpOxciHHlKyhlJaiBLXDc/gFto3rg6qELGflWtiVnsyU5hy1J2Rw8nV+lX1Korwc9mwfRIyaQuOgAOkT6YzHV5ewt2sz6tjMHMexdguHgcgyZhwjP20V43i5UZR5qZHfU2OtRW15HaVg3En5aU+OydwZVVckrKSMj30pGvpUz+VbOFFjJKrRxtrCUrIJSsopKyS4sJbekjOKSEloqaXQgmY5qMp3sSXRSkwlyFFR7/WRHGPvU5ux3NGe/GsN+NYaTaijOrBExGhR8PIz4Wkz4WLRHbaTiBdu5ry/cf/65EYtRZf0vPxMfX3dlfyFVVbHZCqEwE6U4C4qzoTgLpTgbirLAmo9izdUWui7JRSktgNLCc5v2/MImRDM2zOq52pcr7chzkW+RigImT1STJw6TF6rJC4fJC4fRgt3ohd3srT2avCgzelFm9KbM5F3x3GbyokN0J5q27g7U/m/OpZS3RNWEyxKnkJAQjEZjlZqgjIyMKjVG5SIiIqo93mQy0aRJk0sec7FrAlgsFiyWqjU/ZrPZJT/ojtD2KMnrMWcfxdSyn9OvfzW6N2/CqO5RfJN4in/9cJgFD/V17ciaQ9p/BkrsQMwBdbeq/MW+t3aHynNLd/H19pMoCvxzVBfG9a06UMHlWl4LA6fCr29iWjFVG2IfGF33cVxEp2bBvNosmGnDO/KfrSdYuDmF45mFfLMzjW92phEZ4Mmt3aOI7xhB9+jASknUZX+vyqxaX6Zt87Wv296EcvtcTC5osnQ4VI5lFrItOYvNx7PZmpxF8tmq/1m2DPGhd4tgescG07tFEDHB3u4x4iyqs7YN/T+tP8/+ZbBvGUr6LpRT27UpDta/idHDh2s8W2Hx3YMxpg9E9dA6HrtYiIcHIf7edLxwp8OuzZR/5hBkHoSM/XB6D+qZgyj2qn2rHIqJPN9YMv06kOHTjnTvtpyytKJA8cGhqtgdKs1UlSiHikOl4metvC+Pyahc0I9G60tjMRvxPDdtg7YZ8PbQkiNvswlvixEfDxOeZoNT+tkoius+T2rEwwN8rvD3R1W10YNlJdrvZvmjo0zrV+YoO7fZqcik1HPVQxUjDY3aNByK4VzfsAv6iRnMWid3kyeKwQSKgoJzJ5F0RtnX5nyXJU4eHh707NmThISESlMFJCQkcNtt1f933a9fP5YvX15p3+rVq+nVq1fFm+rXrx8JCQmV+jmtXr2a/v37u+BdXBk1pB0kr0fJPKh3KNV6Or4dK3ans/7IWdYezuS6ti5cbLd8GoIO+k96aXeoPPvVTpZsT8VoUJg5pju3dqubKQGqNXgaHPsFUrfB0ke0SRMN7jWKKcDbzKRBLfnDwFh2nMhhyfaTLN+ZRlpuCe+vOcb7a44R7OPB4HahXNe6CUVll7lg7klYPEH7wEeB61+AgU9rf3SdILPAyu7UXHak5LAjJZudJ3LIK6kclKJAhwh/ercIom/LJvRqEUSYXz0YvRTaDkKfhUHPQt4pOPozHPsZjv6MUpRJeOku+HXX+eODWkB4Z2jS+oKtldbf52rK2+GAkhwoyICcFC1Jyk7SHrOOax24y0qqnKaA1gE7vJM2J1dEV4jsiiG0A4FmTwKB1lcelbhSinJuegX36lbizlzaVDd16lQmTJhAr1696NevHx988AEpKSkV8zJNmzaN1NRUPvvsM0AbQTd79mymTp3KpEmT2LhxIx9//HGl0XJPPvkkgwYN4rXXXuO2227j22+/5YcffmDdunWufCu1E9IewG0Tp+hgbyb0a87H647z6vcHGNg65PwoH2fKOaElBSjasHId/T5peueeOEZ0idQ1JoxmuONDeH8QJK/XlmQZ9Iy+MV2Eoij0iAmiR0wQfxvZkR/3Z/D9nnR+OZhBVmEpS7ansmR7KmDio6T19IgJontMIHHRQbQO88XDZIDkjfCfCdo8Op6BcOfH0GbIFcWTX2Ij+WwRR88UsD8tn/1peexLy+NMftXRVZ5mA12bBtI7NojeLYLp0TwIf1eNmqwr/lEQN17bHA5sqYnsXzWPzkElGE7t0NbNy07Stt9TjNp8Qt4h2igzr6ALagjMWg0BKthKtBGOZecei7K1711RplYDcSlGC4S0gZC2ENr+XLLUGQJinJYkC6EXlyZOY8aM4ezZs7z88sukpaXRuXNnVqxYQfPmzQFIS0sjJSWl4vjY2FhWrFjBU089xbvvvktUVBRvv/12xRxOAP3792fRokX89a9/5W9/+xutWrVi8eLFbjWHkxqq9dJQzrhn4gTw+PWt+c/WE+xPy+ObxFTu6OGCeXL2n6s9jOkHfnXXTPd7bpk0lWvSSptd+tvJ8MsMaHk9NOupd1SXZDEZGdElkhFdIrHZHWxLzuanAxn8sO80xzILOXpG2/67TetDZlDgEd9fedr2ISbKyPRpy8be/8ZQHIvfoTP4eprwMhuxO1TKHCpldgc2u0qhtYyzhVZtZFdBKWcLrJzILiYps5CzFxlOrygQG+JD9+hA4mKCiIsOpF2EX8NY4PhiDAaI6MLx0CF0GDECg9kMxTlwagdkHtJqgMq3nBPasjWFZ7TtapbP8wzQEqGg5hDYXHsMaqElTIHN3a72VAhncXnn8MmTJzN58uRqX5s/f36Vfddddx3bt2+/5DXvuusu7rrrLmeE5xJqSDvtSU4ylBbV/eKjNRDk48Hkwa15beUBZnx/gCEdw53/X/jBcyMfOtzi3OvWglsnTeW6j4MjP2hLc3z9EDz6qzavTD1gNhq4pmUTrmnZhGeHtuY/364gpH1vdqfms+NENntPnGWqfR732bRRL9/Z+/Ls2UcoXpkNZF/VvUN8PWjexOfclAD+dIzyp32En/vMgK4nr0Bodb22XaisFIrOarVGhZna8+Ls8zNlO2zacxQwe53fTF7aNSsmbgwFk44TpQqhI/kL4wo+oVhNfljK8rX/+KK66x1RtR4Y0IL/bj3BscxC3lh5kL+P6uy8ixdnQ/K5+braDXfedWvB4VD5y9e73DtpgvNr2Z3cAtnH4Zs/wt2f1csmDV8z3NAulGGdo6DgDOp/X0ZJXo+Kwp72T3Iw6B5uyikhr9hG/gUzX5fYHBdMzqdgMhrw9jDSxMeDJr7n5wZqGuRFiyY+NG/i7bpJShsykwf4R2qbEOKKSOLkIvmeUVgKDmojYdw0cfI0G/nH7Z0Z9+FvLPgtmdt7NKVHjJNGNh35UWsSCG0PwXU71B+0pOlvy3bz1baTGA0Kb49106SpnFeg1ufn05FaE+fP/9DWtquvTu2ARfei5J0EDz+UOz+iS7ub6KJ3XEIIcZXq37+09US+Z1PtyZmrmhHP5fq3CuHOHs1QVXh+yW5sdidNinlopfbY9ibnXK8WVBVe+t9+Fm05gUGBmWO6c3NXN06aysX0hVvf0Z7/+qa2DEc9pOz5L3xyk7Z+WHArmPQjtKv7nwMhhHAFSZxc5Hzi5L4dxMu9cHMHgrzNHEjP55N1x6/+gvYyOHxuJtc6TpxUVWVJkoGFm7V5mt4c3U3fKQdqq9tYbXg+wPIntJFo9YWjjE4nF2L69o/aSKw28TDpJ20YvRBCNBCSOLlIReKU4d41TgDBPh48P6IDALN+OMyJrJpPPV+tE5u0eV68giG6z9UHWEOqqvLqykOsTdd+rF+7syu3x7lgtKCrXf9Xbd4reyksHq/NjePuCs5g/HI0rc+cq2kc+Azcs0hrghRCiAZEEicXqUicspPAVqxrLDVxV89m9I0Npthm5/++3VOx1M0VOfi99tgmvs6GJKuqyj/+t59PNiQD8I/bOjK6l/vMxF0rBgPc/j5EdtdGPS0coy2d4K6S1sHcazEkraXMYKHsjk/gxr/JcHQhRIMkiZOLWE3+qF7BgKqNrHNziqLwz9u7YDYq/HzwDN8mnrryix1apT3WUb8WVVV5afk+Pj7XzDi6pZ0xvephTdOFPLy1Ghu/KG3ZinnDIS9N76gqczhg7b/g01ugIB01pB1r276I6gazxAshhKtI4uQqioIacm650nrQzwmgdZgvj1/fBoDnl+7mSEZ+7S9y9qg2a7HBBK1ucHKEVamqyovL9jJ/QxKg1TQNCL+K2jJ34h8J930DfpFw5gDMu6n6maD1UJgJC++Gn/4OqgO63UPZA6vJ96rnCasQQlyGJE4uVDER5pkD+gZSC49d34p+LZtQVGrnjwu2U1R6ucXHfqe8ma75AG1mYRdyOFT++s0ePtuYjKLA63d2rf81Tb8X2g4eXKnNyJydpI1Wy9D552nvNzDnGm3STpMn3DobRr0HHj76xiWEEHVAEidXCtXWrNP9g64WTEYD/76nO2F+Fg5nFPD8kt216+9UPg2Biye9LLM7mLZkN1/8loKiwBt3dWN073rap+lyglrAAyshtAPkp2nNdqd21H0cBWfgP/fBf+/XlusIba+NmusxQZvEUwghGgFJnFyoPtY4AYT5efLOPXEYDQrfJJ5i4eaUy58E2vpY5bOFtx3msvhKbHYeXbCdxVu1eZreGt2Nu3o2sJqm3/OPhAdWQFQcFGfBJ8Nh03taPyNXU1XY/RW82wf2fastEjvoWXhkrbZ4qxBCNCKSOLmQWj5/TfZxbaXxeqRvyyY8O0yL/6Vl+9h9MvfyJx35QZstPKQdBLd0SVw5RaWM/+g3fth/Gg+TgTnje9bPKQeuhHcw3LcMWt0IZcWw8jltpvGsY667Z/IGmDdCW0OvOAvCu8DDP8MNfwWTxXX3FUIINyWJkyv5hIFnoNZ59uxhvaOptYcHtmRIh3BK7Q4mL9xGRv5lkr+KZjrXjKY7lVPM3XM3si05G39PEwse6stNnSNcci+35ekP934NN78FZh9IXg/vDYDNHzq39il1G3x+u9YsmLIBjBa4/gUtaYrs5rz7CCFEPSOJkyspCoRpE0vWp35O5QwGhTfv7kZ0sBcnsooZ+8EmTuddJHmqNFu48/s37U/L4873NnA4o4Bwfwv/fbQ/fWKDnX6fekFRoPdDMHkDtBgItiJY8QzMHQC/faA1mV4JWzHsWwYLx8KHN8DRn7TRkb0ehCcT4bo/g1EW1hVCNG6SOLlaaP3s51QuwNvMgof6EhXgybEzhYz9YBNpudVM6Hnit3OzhQdBs95OjeG/W09w+5z1pOWW0CrUh6//2J92EX5OvUe9FNRCa7ob/oZW+5SxD75/Ft5sD0v/qC3Xcrkm4tJC2P8dfPUQvNEa/jMBDn0PigG6jYPHt8LImeBfj5atEUIIFzLpHUCDF3quxqmeJk4AzZv4sPiRftzz4SaOZxYy5v1NLJzUl2ZB3ucPOnxu0svWQ8HonB+r4lJtFvP/bjsJwMA2Ibw9No4gHw+nXL9BMBig78PQ9W7Y9V/YNk9LoHYu1DYA3wgtyQpqrk0RkXcKck9A7kltZvILBURDp1EQdx+Etq3rdyOEEG5PEidXq+c1TuWig7215OmDTaRkFTHm/U0sevgaooPPJU+HVmuPThpNd/RMAY99sZ0D6fkYFHhqSFseu741BoMMe6+WV5CWQPWZBCe3wrb5sH85WHOhIF3bTmyq/lz/Zlqy1Ol2aNpTphYQQohLkMTJ1crncso6pjWbmD31jecqNA30YvEj1zDuw984nlnIqHfX8+KtnbglxoZyZr/WvHOVs4WX2OzM35DEOz8eprDUToivhbfHdqd/6xAnvYsGTlEgure23TYbirO1iTOzkyAnGUpywb+pVrMU0EzbPAMkWRJCiBqSxMnV/CK0D6aSXDh7BCI66x3RVYkM8GLRw9dw/yebOZCezxNf7iA7ahP3AzTrow2ZvwKqqvLdrjReW3mAk9laH6q+scG8c08cYf71N9nUlaJo3w/vYGjaQ+9ohBCiQZDO4a6mKOdrnep5c125cH9Plj1+LU8PbYuH0UDTM78CsNXSG7ujduvEORwqG4+e5Y73NvCnL3dwMruYCH9P3ry7G19OukaSJiGEEG7FpYlTdnY2EyZMICAggICAACZMmEBOTs4lz1FVlenTpxMVFYWXlxeDBw9m7969lY4ZPHgwiqJU2saOHevCd3KVGljiBOBhMvCnG9vw/WO9GWjcB8Df9kbR558/8PR/drJidxr5JbZqz7WW2fnlYAbPL91N3xk/cs+Hm9iRkoOX2cjUoW35+ZnB3NmzmfRnEkII4XZc2lQ3btw4Tp48ycqV2sSIDz/8MBMmTGD58uUXPef111/nrbfeYv78+bRt25Z//OMfDB06lIMHD+Lnd34I+qRJk3j55Zcrvvby8nLdG7laFWvW7dc3DhdoVZgIWCmyhJGmtCKnsJSvt5/k6+0nMRsVukcHYjEZKbU7sJ3bkjKLKLCeXzzYz2JiZLcopgxpQ7jUMAkhhHBjLkuc9u/fz8qVK9m0aRN9+/YF4MMPP6Rfv34cPHiQdu3aVTlHVVVmzZrFCy+8wB133AHAp59+Snh4OAsXLuSRRx6pONbb25uIiJrNGm21WrFarRVf5+XlAWCz2bDZqq8VuRrl1yx/VEI7YALUtJ2UueB+ejIcXIkR8OwwjA03DWZ7Sg4/HTjDzwfPcPxsEVuSsqs9L8zPwo3tQxnaIYy+scF4mLTKz6v9fvy+7EXdkbLXj5S9fqTs9ePMsq/NNRRVVWvXKaWGPvnkE6ZOnVqlaS4wMJCZM2fywAMPVDnn2LFjtGrViu3btxMXF1ex/7bbbiMwMJBPP/0UoKL5TlVVwsPDGT58OC+++GKlGqkLTZ8+nZdeeqnK/oULF+Lt7V3NGc5lKivk5t1/BGBFlznYTL4uv2edUFWG7HsWn9IMfot9kvTAnpVeziiG5AIFBTAZwKhom59ZpakPSEucEEIId1BUVMS4cePIzc3F39//kse6rMYpPT2dsLCwKvvDwsJIT0+/6DkA4eHhlfaHh4eTnJxc8fX48eOJjY0lIiKCPXv2MG3aNHbu3ElCQkK11502bRpTp06t+DovL4/o6Gji4+MvW0BXwmazkZCQwNChQzGbtSUq1BOvoeQkEd8lAjV2kNPvqYuzhzEnZqAazPS46ynw0D8hrK7sRd2QstePlL1+pOz148yyL2+JqolaJ04Xq7250JYtWwBQqpkbRlXVavdf6Pev//6cSZMmVTzv3Lkzbdq0oVevXmzfvp0ePaoOu7ZYLFgsVVdyN5vNLv1Br3T9yK6Qk4TpzF5oe6PL7lmnjv8MgNJiAGafIJ2DqczV31txcVL2+pGy14+UvX6cUfa1Ob/WidPjjz9+2RFsLVq0YNeuXZw+fbrKa2fOnKlSo1SuvM9Seno6kZGRFfszMjIueg5Ajx49MJvNHD58uNrEyS1EdoP9yyBtp96ROM+hc8ustHHObOFCCCGEu6t14hQSEkJIyOVnce7Xrx+5ubls3ryZPn36APDbb7+Rm5tL//79qz2nvPktISGhoo9TaWkpa9as4bXXXrvovfbu3YvNZquUbLmdyO7aY/ouXcNwGms+JG/QnreJ1zcWIYQQoo64bB6nDh06cNNNNzFp0iQ2bdrEpk2bmDRpEiNHjqw0oq59+/YsXboU0JropkyZwiuvvMLSpUvZs2cPEydOxNvbm3HjxgFw9OhRXn75ZbZu3UpSUhIrVqzg7rvvJi4ujgEDBrjq7Vy9yK7aY+ZhsBboG4szHFsDDhsExUKTVnpHI4QQQtQJl87j9MUXX/DEE08QH6/VSNx6663Mnj270jEHDx4kNze34us///nPFBcXM3nyZLKzs+nbty+rV6+uGDHn4eHBjz/+yL///W8KCgqIjo7m5ptv5sUXX8RoNLry7Vwd3zDwi4T8NDi9F2L66h3R1Tl8rpmu7TBZ50wIIUSj4dLEKTg4mAULFlzymN/PhqAoCtOnT2f69OnVHh8dHc2aNWucFWLdiuiqJU5pO+t34qSqcPjcCMY2Q/WNRQghhKhDslZdXYrspj2m1/MO4um7tQTQ7A3Nr9U7GiGEEKLOSOJUl8r7OdX3kXXlo+laDgazLJEihBCi8ZDEqS6V1zhlHIAy66WPdWeHtLUHaSvTEAghhGhcJHGqSwHR4BmojUarrwv+FpyB1G3ac5mGQAghRCMjiVNdUpQL+jnV0/mcjiQAqtbR3T9K72iEEEKIOiWJU12r7/2cKprpbtI3DiGEEEIHkjjVtfIZxNPqYY1TWSkc+Ul7LomTEEKIRkgSp7oWca7G6fQecNj1jaW2UjZCaT74hEJUnN7RCCGEEHVOEqe61qQVmH3AVgRnj+gdTe1cuKivQX50hBBCND7y6VfXDEaI6Kw9r2/9nCr6N8loOiGEEI2TJE56KB9ZV58Sp8wjkHUUDGZoeb3e0QghhBC6kMRJDxH1cGRdeW1TiwHg6a9vLEIIIYROJHHSw4VzOf1ukWO3dfhc/yYZTSeEEKIRk8RJD6HttSavklzISdY7mssryYXkDdpzmS1cCCFEIyaJkx5MHhDeUXteH5rrjv4EjjJo0kYbFSiEEEI0UpI46aV8HqSTW/SNoybKpyGQRX2FEEI0cpI46SWmv/aYvFHfOC7HXgaHV2vPJXESQgjRyEnipJfm/bTHtESwFugayiWlbISis+AVdD7ZE0IIIRopSZz0EhgDAdFa3yF3bq7bv1x7bDcCjCZ9YxFCCCF05tLEKTs7mwkTJhAQEEBAQAATJkwgJyfnkucsWbKEYcOGERISgqIoJCYmVjnGarXypz/9iZCQEHx8fLj11ls5efKka96EK8Wcq3VKcdPmOlWFA99pzzvcom8sQgghhBtwaeI0btw4EhMTWblyJStXriQxMZEJEyZc8pzCwkIGDBjAq6++etFjpkyZwtKlS1m0aBHr1q2joKCAkSNHYrfXs0Vzm5f3c9qgbxwXc2o75KVqa+vJbOFCCCEELmt72b9/PytXrmTTpk307dsXgA8//JB+/fpx8OBB2rVrV+155YlVUlJSta/n5uby8ccf8/nnnzNkyBAAFixYQHR0ND/88APDhtWjDszNB2iPJ7dAmRVMFn3j+b3yZro2Q8HsqW8sQgghhBtwWeK0ceNGAgICKpImgGuuuYaAgAA2bNhw0cTpcrZt24bNZiM+/vxEjFFRUXTu3JkNGzZUmzhZrVasVmvF13l5eQDYbDZsNtsVxXEp5de87LUDWmDyDkEpyqTsxFbUZn2cHssVU1VM+5ahAGVtR6C6oJxcocZlL5xOyl4/Uvb6kbLXjzPLvjbXcFnilJ6eTlhYWJX9YWFhpKenX9V1PTw8CAoKqrQ/PDz8otedMWMGL730UpX9q1evxtvb+4pjuZyEhITLHtPb3IIoMjm0eh6HIzJdFktt+RWnckPWUeyKiVXHVMqSV+gdUq3UpOyFa0jZ60fKXj9S9vpxRtkXFRXV+NhaJ07Tp0+vNgm50JYt2igxRVGqvKaqarX7r9alrjtt2jSmTp1a8XVeXh7R0dHEx8fj7+/8BWttNhsJCQkMHToUs9l8yWMNm1MgYSvtvbNpM2KE02O5UoZ1b8IBUFoOJv6WO/UOp8ZqU/bCuaTs9SNlrx8pe/04s+zLW6JqotaJ0+OPP87YsWMveUyLFi3YtWsXp0+frvLamTNnCA8Pr+1tK0RERFBaWkp2dnalWqeMjAz6969+niGLxYLFUrX/kNlsdukPeo2uHzsQAMPJzRiMBjAYXRZPrRz8HwCGTrdhqId/DFz9vRUXJ2WvHyl7/UjZ68cZZV+b82udOIWEhBASEnLZ4/r160dubi6bN2+mTx+t785vv/1Gbm7uRROcmujZsydms5mEhARGjx4NQFpaGnv27OH111+/4uvqJqILePiBNQ9O74HIbnpHBNnJkL4LFIM2f5MQQgghABdOR9ChQwduuukmJk2axKZNm9i0aROTJk1i5MiRlTqGt2/fnqVLl1Z8nZWVRWJiIvv27QPg4MGDJCYmVvRfCggI4KGHHuLpp5/mxx9/ZMeOHdx777106dKlYpRdvWIwQsy5DvTusvxK+dxNMf3B5/JJshBCCNFYuHQepy+++IIuXboQHx9PfHw8Xbt25fPPP690zMGDB8nNza34etmyZcTFxXHzzTcDMHbsWOLi4pg7d27FMTNnzmTUqFGMHj2aAQMG4O3tzfLlyzEa3aSZq7Yq5nNar28c5cqnIZBJL4UQQohKXLqGRnBwMAsWLLjkMaqqVvp64sSJTJw48ZLneHp68s477/DOO+9cbYjuoXw+p+QN2mzdLug8X2MFGZCySXveYaR+cQghhBBuSNaqcwdRcWC0QFEmnD2ibywH/geoENUDAprpG4sQQgjhZiRxcgcmCzTrrT3Xu7luz9fao9Q2CSGEEFVI4uQump9b8FfPdetyUiDpV+15l7v1i0MIIYRwU5I4uYuKDuI6jqzbuVh7bDEQAmP0i0MIIYRwU5I4uYtmfUAxQm4KZCfV/f1VFXZ+qT3vPq7u7y+EEELUA5I4uQuLL0Sfm8/pwP/q/v4nNkPWUTB7Q4db6/7+QgghRD0giZM76Xib9rjv27q/986F2mOHW7UkTgghhBBVSOLkTjqeq+k58Rvknaq7+9qKYc+52du731N39xVCCCHqGUmc3Il/1Pnmun3L6u6+B1eANRf8m0GLQXV3XyGEEKKekcTJ3XQcpT3WZXPdzkXaY7cxYJAfCSGEEOJi5FPS3ZQ316VshPx0198v/zQc+VF73k2a6YQQQohLkcTJ3QQ0OzeLuHp+sV1X2v0fUO3aPUPauP5+QgghRD0miZM7Kh9dt/cb195HVSHx3NxNUtskhBBCXJYkTu6oPHFKXg8FGa67T/ouyNirLTDc+Q7X3UcIIYRoICRxckeBMRDVA625zoWj6zbM1h7b3wxeQa67jxBCCNFASOLkrjqN0h5dNbou6xjs+Up7PuAJ19xDCCGEaGAkcXJX5c11SeugMNP51183E1QHtB4KUXHOv74QQgjRAEni5K6CWkBkdy25cfboutyT5zuFD3rGudcWQgghGjBJnNxZRXPdN8697vq3wWGDFgMh5hrnXlsIIYRowCRxcmcdRwEKHPsFTiU655oFGbD9U+35wKedc00hhBCikXBp4pSdnc2ECRMICAggICCACRMmkJOTc8lzlixZwrBhwwgJCUFRFBITE6scM3jwYBRFqbSNHTvWNW9CT8Gx0OUu7fmPLznnmhtnQ1kJNO0FLQc755pCCCFEI+HSxGncuHEkJiaycuVKVq5cSWJiIhMmTLjkOYWFhQwYMIBXX331ksdNmjSJtLS0iu399993Zuju4/oXwGCGoz/BsTVXd62iLNjysfZ80DOgKFcfnxBCCNGImFx14f3797Ny5Uo2bdpE3759Afjwww/p168fBw8epF27dtWeV55YJSUlXfL63t7eRERE1CgWq9WK1Wqt+DovLw8Am82GzWar0TVqo/yaTrm2XzMMPe7HuPUjHD9Mxz5x1RUnPIaNczCWFqCGdaYs9kZwwXvXm1PLXtSKlL1+pOz1I2WvH2eWfW2uoaiqql71HavxySefMHXq1CpNc4GBgcycOZMHHnjgkucnJSURGxvLjh076N69e6XXBg8ezN69e1FVlfDwcIYPH86LL76In59ftdeaPn06L71Utalr4cKFeHt71+p96cFiy2XIvmcwOaxsjv0TaYG9r+gaN+x/Dg97IVtaPMapoL4uiFQIIYSof4qKihg3bhy5ubn4+/tf8liX1Tilp6cTFhZWZX9YWBjp6elXde3x48cTGxtLREQEe/bsYdq0aezcuZOEhIRqj582bRpTp06t+DovL4/o6Gji4+MvW0BXwmazkZCQwNChQzGbzU65phKYBOvepHfe95SNfQEMtfjWOcowLrwTg70QNawj3e/5P7objE6Jy924ouxFzUjZ60fKXj9S9vpxZtmXt0TVRK0Tp4vV3lxoy5YtACjVNCmpqlrt/tqYNGlSxfPOnTvTpk0bevXqxfbt2+nRo0eV4y0WCxaLpcp+s9ns0h90p17/2imwbR7K2SOY9/wHet5f83NXTdfWvfPwRbn7U8wWT+fE5MZc/b0VFydlrx8pe/1I2evHGWVfm/NrnTg9/vjjlx3B1qJFC3bt2sXp06ervHbmzBnCw8Nre9tL6tGjB2azmcOHD1ebODUInv5ah+5Vz8Mvr0LX0WD2uvx5e5dqI+kARs2B0LaujVMIIYRowGqdOIWEhBASEnLZ4/r160dubi6bN2+mT58+APz222/k5ubSv3//2kd6CXv37sVmsxEZGenU67qdXg/Bpvcg9wT88BIM+ydcqskt4wB885j2fMCT55dxEUIIIcQVcdl0BB06dOCmm25i0qRJbNq0iU2bNjFp0iRGjhxZaURd+/btWbp0acXXWVlZJCYmsm/fPgAOHjxIYmJiRb+oo0eP8vLLL7N161aSkpJYsWIFd999N3FxcQwYMMBVb8c9mD3hxv/Tnv/2Hnx+O+RXrdUDoCQPFo8HW6E2Q/gN/1d3cQohhBANlEvncfriiy/o0qUL8fHxxMfH07VrVz7//PNKxxw8eJDc3NyKr5ctW0ZcXBw333wzAGPHjiUuLo65c+cC4OHhwY8//siwYcNo164dTzzxBPHx8fzwww8YjQ2zw3MlXUfDbXPA7A3H18DcAXDkh/Ov56bCprkwfwScPQL+TeGueWB02TgAIYQQotFw6adpcHAwCxYsuOQxv58NYeLEiUycOPGix0dHR7NmzVVOBFnfxY2HZr3gvw9Axl5YcCd0G6clSic3nz/O7A2jPwPfUP1iFUIIIRoQqYaor0LbwaQfYdULsPVj2Lnw/GvR12j9mTqNAv8o3UIUQgghGhpJnOozsxeMfAtaXQ+7FkOLQdBhpCRLQgghhItI4tQQdLhF24QQQgjhUi7tHC6EEEII0ZBI4iSEEEIIUUOSOAkhhBBC1JAkTkIIIYQQNSSJkxBCCCFEDTXKUXXlk27m5eW55Po2m42ioiLy8vJktew6JmWvHyl7/UjZ60fKXj/OLPvyfOD3k3JXp1EmTvn5+YA2C7kQQgghBGj5QUBAwCWPUdSapFcNjMPh4NSpU/j5+aEoitOvn5eXR3R0NCdOnMDf39/p1xcXJ2WvHyl7/UjZ60fKXj/OLHtVVcnPzycqKgqD4dK9mBpljZPBYKBZs2Yuv4+/v7/8IulEyl4/Uvb6kbLXj5S9fpxV9peraSonncOFEEIIIWpIEichhBBCiBqSxMkFLBYLL774IhaLRe9QGh0pe/1I2etHyl4/Uvb60avsG2XncCGEEEKIKyE1TkIIIYQQNSSJkxBCCCFEDUniJIQQQghRQ5I4CSGEEELUkCROTjZnzhxiY2Px9PSkZ8+e/Prrr3qH1CjMmDGD3r174+fnR1hYGKNGjeLgwYN6h9XozJgxA0VRmDJlit6hNBqpqance++9NGnSBG9vb7p37862bdv0DqvBKysr469//SuxsbF4eXnRsmVLXn75ZRwOh96hNThr167llltuISoqCkVR+Oabbyq9rqoq06dPJyoqCi8vLwYPHszevXtdFo8kTk60ePFipkyZwgsvvMCOHTsYOHAgw4cPJyUlRe/QGrw1a9bw2GOPsWnTJhISEigrKyM+Pp7CwkK9Q2s0tmzZwgcffEDXrl31DqXRyM7OZsCAAZjNZr7//nv27dvHm2++SWBgoN6hNXivvfYac+fOZfbs2ezfv5/XX3+dN954g3feeUfv0BqcwsJCunXrxuzZs6t9/fXXX+ett95i9uzZbNmyhYiICIYOHVqxLq3TqcJp+vTpoz766KOV9rVv31597rnndIqo8crIyFABdc2aNXqH0ijk5+erbdq0URMSEtTrrrtOffLJJ/UOqVH4y1/+ol577bV6h9Eo3XzzzeqDDz5Yad8dd9yh3nvvvTpF1DgA6tKlSyu+djgcakREhPrqq69W7CspKVEDAgLUuXPnuiQGqXFyktLSUrZt20Z8fHyl/fHx8WzYsEGnqBqv3NxcAIKDg3WOpHF47LHHuPnmmxkyZIjeoTQqy5Yto1evXtx9992EhYURFxfHhx9+qHdYjcK1117Ljz/+yKFDhwDYuXMn69atY8SIETpH1rgcP36c9PT0Sp+9FouF6667zmWfvY1ykV9XyMzMxG63Ex4eXml/eHg46enpOkXVOKmqytSpU7n22mvp3Lmz3uE0eIsWLWL79u1s2bJF71AanWPHjvHee+8xdepUnn/+eTZv3swTTzyBxWLhvvvu0zu8Bu0vf/kLubm5tG/fHqPRiN1u55///Cf33HOP3qE1KuWfr9V99iYnJ7vknpI4OZmiKJW+VlW1yj7hWo8//ji7du1i3bp1eofS4J04cYInn3yS1atX4+npqXc4jY7D4aBXr1688sorAMTFxbF3717ee+89SZxcbPHixSxYsICFCxfSqVMnEhMTmTJlClFRUdx///16h9fo1OVnryROThISEoLRaKxSu5SRkVElExau86c//Ylly5axdu1amjVrpnc4Dd62bdvIyMigZ8+eFfvsdjtr165l9uzZWK1WjEajjhE2bJGRkXTs2LHSvg4dOvD111/rFFHj8eyzz/Lcc88xduxYALp06UJycjIzZsyQxKkORUREAFrNU2RkZMV+V372Sh8nJ/Hw8KBnz54kJCRU2p+QkED//v11iqrxUFWVxx9/nCVLlvDTTz8RGxurd0iNwo033sju3btJTEys2Hr16sX48eNJTEyUpMnFBgwYUGXajUOHDtG8eXOdImo8ioqKMBgqf4QajUaZjqCOxcbGEhERUemzt7S0lDVr1rjss1dqnJxo6tSpTJgwgV69etGvXz8++OADUlJSePTRR/UOrcF77LHHWLhwId9++y1+fn4VNX8BAQF4eXnpHF3D5efnV6UfmY+PD02aNJH+ZXXgqaeeon///rzyyiuMHj2azZs388EHH/DBBx/oHVqDd8stt/DPf/6TmJgYOnXqxI4dO3jrrbd48MEH9Q6twSkoKODIkSMVXx8/fpzExESCg4OJiYlhypQpvPLKK7Rp04Y2bdrwyiuv4O3tzbhx41wTkEvG6jVi7777rtq8eXPVw8ND7dGjhwyHryNAtdu8efP0Dq3RkekI6tby5cvVzp07qxaLRW3fvr36wQcf6B1So5CXl6c++eSTakxMjOrp6am2bNlSfeGFF1Sr1ap3aA3Ozz//XO3f9/vvv19VVW1KghdffFGNiIhQLRaLOmjQIHX37t0ui0dRVVV1TUomhBBCCNGwSB8nIYQQQogaksRJCCGEEKKGJHESQgghhKghSZyEEEIIIWpIEichhBBCiBqSxEkIIYQQooYkcRJCCCGEqCFJnIQQQgghakgSJyGEEEKIGpLESQjRaA0ePJgpU6boHYYQoh6RxEkIIYQQooZkrTohRKM0ceJEPv3000r7jh8/TosWLfQJSAhRL0jiJIRolHJzcxk+fDidO3fm5ZdfBiA0NBSj0ahzZEIId2bSOwAhhNBDQEAAHh4eeHt7ExERoXc4Qoh6Qvo4CSGEEELUkCROQgghhBA1JImTEKLR8vDwwG636x2GEKIekcRJCNFotWjRgt9++42kpCQyMzNxOBx6hySEcHOSOAkhGq1nnnkGo9FIx44dCQ0NJSUlRe+QhBBuTqYjEEIIIYSoIalxEkIIIYSoIUmchBBCCCFqSBInIYQQQogaksRJCCGEEKKGJHESQgghhKghSZyEEEIIIWpIEichhBBCiBqSxEkIIYQQooYkcRJCCCGEqCFJnIQQQgghakgSJyGEEEKIGvp/SBwicpxozJcAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 600x250 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(6,2.5))\n",
    "for i in range(y.shape[1]):\n",
    "    plt.plot(t, y[:,i], label=f'y[{i}]')\n",
    "plt.xlabel('t')\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.savefig('images/pendulum-time-plot.png')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate dydt values\n",
    "dydt = pendulum_dydt(t, y.T).T\n",
    "assert dydt.shape == (100, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dynopt.models.models import NonLinearModel\n",
    "\n",
    "# Labels for inputs and outputs\n",
    "x_names = ['theta', 'theta_dot']\n",
    "y_names = ['theta_dot', 'theta_ddot']\n",
    "\n",
    "# Choose input features including non-linear terms\n",
    "x_features = ['x1', 'sin(x0)']\n",
    "\n",
    "# Initialize model\n",
    "model = NonLinearModel(x_names, y_names, x_features=x_features)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      theta  theta_dot\n",
      "0  0.100000   0.000000\n",
      "1  0.097595  -0.047109\n",
      "2  0.090802  -0.087513\n",
      "3  0.080359  -0.119943\n",
      "4  0.067106  -0.143621\n"
     ]
    }
   ],
   "source": [
    "# Prepare pandas dataframes for inputs and outputs\n",
    "y = pd.DataFrame(y, columns=x_names)\n",
    "dydt = pd.DataFrame(dydt, columns=y_names)\n",
    "print(y.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   theta_dot  theta_ddot\n",
      "0   0.000000   -0.499167\n",
      "1  -0.047109   -0.440093\n",
      "2  -0.087513   -0.365874\n",
      "3  -0.119943   -0.281420\n",
      "4  -0.143621   -0.191657\n"
     ]
    }
   ],
   "source": [
    "print(dydt.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     x1  sin(x0)\n",
      "y0  1.0     -0.0\n",
      "y1 -1.0     -5.0\n"
     ]
    }
   ],
   "source": [
    "# Fit the linear model\n",
    "model.fit(y, dydt)\n",
    "print(model.coef_.round(5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "y0   -0.0\n",
      "y1    0.0\n",
      "dtype: float64\n"
     ]
    }
   ],
   "source": [
    "print(model.intercept_.round(5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dynopt.models.models import SparseNonLinearModel\n",
    "\n",
    "# Labels for inputs and outputs\n",
    "x_names = ['theta', 'theta_dot']\n",
    "y_names = ['theta_dot', 'theta_ddot']\n",
    "\n",
    "# Choose additional non-linear features for model identification\n",
    "custom_features = ['sin(x0)', 'cos(x0)', 'sin(x0)**2', \n",
    "                   'cos(x0)**2', 'sin(x0)*cos(x0)']\n",
    "\n",
    "# Initialize SINDy model\n",
    "model = SparseNonLinearModel(x_names, y_names, \n",
    "                             custom_features=custom_features,\n",
    "                             poly_order=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prepare pandas dataframes for inputs and outputs\n",
    "y = pd.DataFrame(y, columns=x_names)\n",
    "dydt = pd.DataFrame(dydt, columns=y_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3\n"
     ]
    }
   ],
   "source": [
    "# Sparse identification of non-linear dynamics (SINDy)\n",
    "threshold = 0.2  # Sparsity parameter\n",
    "model.fit(y, dydt, threshold=threshold)\n",
    "print(model.n_params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     x1  sin(x0)\n",
      "y0  1.0      0.0\n",
      "y1 -1.0     -5.0\n"
     ]
    }
   ],
   "source": [
    "print(model.coef_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "torch",
   "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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}