example_use.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Training symboloc regression models\n",
    "This notebook illustrates how to use  scripts developed to train symbolic regression models. For more information refer to the [main manuscript](https://doi.org/10.1039/D2DD00027J)  \n",
    "\n",
    "The following examples were run with the following dependencies:\n",
    "* numpy == 1.20.0\n",
    "* pandas == 1.2.1\n",
    "* sympy == 1.9\n",
    "* scikit-learn == 1.0.1\n",
    "* autofeat == 2.0.10  \n",
    "\n",
    "If you have unexpected errors you might try to reproduce this environment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import workflows as wf"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training data\n",
    "\n",
    "For this example we generate artifical training data from a arbitrarily chosen expression. In your application, however, you will replace this function for a function opening an existing dataset from, e.g. a pickle file, or a csv. Your training data should be in a pandas dataframe of shape (n samples) x (m predictors + target)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# -------------- TRAINING DATA -------------------\n",
    "\n",
    "def generate_training_data():\n",
    "    \"\"\"Example function to generate toy training set\"\"\"\n",
    "\n",
    "    npoints = 50\n",
    "    t_coords = np.linspace(243, 330, npoints) #temperature samples\n",
    "    c_coords = np.linspace(0.1, 3, npoints) #concentration samples\n",
    "    b1, b2 = 1, -2e-3 #coefficients\n",
    "\n",
    "    #training dataframe\n",
    "    train_data = pd.DataFrame(data = np.array([np.repeat(t_coords,npoints**2), np.hstack((npoints**2)*(c_coords,))]).T, \n",
    "                            columns=['T', 'c'])\n",
    "\n",
    "    #apply expression to be found via symbolic regression: k = b1*c*t**1.5 + b2*t*c**1.5\n",
    "    train_data['k'] = b1*train_data['c']*train_data['T'].pow(0.5) + b2*train_data['T']*train_data['c'].pow(1.5)\n",
    "    \n",
    "    return train_data\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting hyperparameters\n",
    "The workflow can be conceptually divided into:\n",
    "* A feature generation step, where hundreds of features are generated from the initial predictors by combining these using non-linear transformations.\n",
    "* A feature selection step, where  most candidate features are filtered out due to multiple criteria, e.g. close-to-zero coefrficients, non-sensical dimensions, poor correlation to target, high-correlation to simpler features, etc.\n",
    "\n",
    "More details about the steps consult the [Autofeat documentation](https://github.com/cod3licious/autofeat).\n",
    "\n",
    "We first set all hyperparamters, for the feature generation and selection steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# --------------  HYPER-PARAMETERS FEATURE GENERATION -------------------\n",
    "\n",
    "#how to scale data. Supported 'standard_nomean', 'standard', 'none'\n",
    "SCALING_TYPE = 'standard_nomean' \n",
    "\n",
    "#whether to leave intercept to vary freely (True) or constrain its value to y0 = 0 (False).\n",
    "FIT_INTERCEPT = False\n",
    "\n",
    "# Autofeat hyperparameter.\n",
    "# number of times features are combined to obtain ever more complex features.\n",
    "# example FEATENG_STEPS = 3 with sqrt transformations will find terms like sqrt(sqrt(sqrt(x)))\n",
    "FEATENG_STEPS = 5\n",
    "\n",
    "# Autofeat hyperparameter.\n",
    "# Units of predictors. Keys must match column names in dataframe. \n",
    "# Ignored predictors are assumed to be dimensionless.\n",
    "UNITS = {\"T\": \"1/K\",\n",
    "        \"c\": \"mol/kg\"}\n",
    "\n",
    "# Autofeat hyperparameter.\n",
    "# Number of iterations for filtering out generated features.\n",
    "FEATSEL_RUNS = 3\n",
    "\n",
    "# Autofeat hyperparameter.\n",
    "# Set of non-linear transformations to be applied to initial predictors.\n",
    "# Autofeat documents these possible transformations: [\"1/\", \"exp\", \"log\", \"abs\", \"sqrt\", \"^2\", \"^3\", \"1+\", \"1-\", \"sin\", \"cos\", \"exp-\", \"2^\"].\n",
    "# Autofeat throws an error when using a single transformation. \n",
    "# Repeat your transformation as a workaround if you only want o use one.\n",
    "TRANSFORMATIONS = [\"sqrt\", \"sqrt\"]\n",
    "\n",
    "\n",
    "# --------------  HYPER-PARAMETERS FEATURE SELECTION -------------------\n",
    "\n",
    "# n-standard deviations criterion to choose optimal alpha from Cross Validation. \n",
    "# Higher STD_ALPHA lead to sparser solutions.\n",
    "STD_ALPHA = 1 \n",
    "\n",
    "#t-statistic rejection threshold. Coefficients with t-statistic < REJECTION_THR are rejected.\n",
    "REJECTION_THR = 2 \n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run workflows\n",
    "Workflows are essentially pre-defined recipes of training steps:\n",
    "* standardization\n",
    "* feature generation\n",
    "* determination of coecfficients and t-values\n",
    "* feature selection\n",
    "* final training\n",
    "* generate trained model objects\n",
    "\n",
    "Below we instantiate a workflow object, with the required parameters and data. Then we run the workflow method that execute the steps in the recipe. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[AutoFeat] Warning: This just calls fit_transform() but does not return the transformed dataframe.\n",
      "[AutoFeat] It is much more efficient to call fit_transform() instead of fit() and transform()!\n",
      "[AutoFeat] Applying the Pi Theorem\n",
      "[AutoFeat] The 5 step feature engineering process could generate up to 1354 features.\n",
      "[AutoFeat] With 125000 data points this new feature matrix would use about 0.68 gb of space.\n",
      "[feateng] Step 1: transformation of original features\n",
      "[feateng] Generated 2 transformed features from 2 original features - done.\n",
      "[feateng] Step 2: first combination of features\n",
      "[feateng] Generated 6 feature combinations from 6 original feature tuples - done.\n",
      "[feateng] Step 3: transformation of new features\n",
      "[feateng] Generated 6 transformed features from 6 original features - done.\n",
      "[feateng] Step 4: combining old and new features\n",
      "[feateng] Generated 36 feature combinations from 48 original feature tuples - done.\n",
      "[feateng] Step 5: combining new features\n",
      "[feateng] Generated 58 feature combinations from 66 original feature tuples - done.\n",
      "[feateng] Generated altogether 108 new features in 5 steps\n",
      "[feateng] Removing correlated features, as well as additions at the highest level\n",
      "[feateng] Generated a total of 8 additional features\n",
      "[featsel] Scaling data...done.\n",
      "[featsel] Feature selection run 1/3\n",
      "[featsel] Feature selection run 2/3\n",
      "[featsel] Feature selection run 3/3\n",
      "[featsel] 5 features after 3 feature selection runs\n",
      "[featsel] 3 features after correlation filtering\n",
      "[featsel] 2 features after noise filtering\n",
      "[AutoFeat] Computing 1 new features.\n",
      "[AutoFeat]     1/    1 new features ...done.\n",
      "[AutoFeat] Final dataframe with 3 feature columns (1 new).\n",
      "[AutoFeat] Training final regression model.\n",
      "[AutoFeat] Trained model: largest coefficients:\n",
      "-2.127543409731164\n",
      "13.618914 * c\n",
      "0.000061 * T**(3/2)*sqrt(T*c)\n",
      "[AutoFeat] Final score: 0.9994\n",
      "[AutoFeat] Applying the Pi Theorem\n",
      "[AutoFeat] Computing 1 new features.\n",
      "[AutoFeat]     1/    1 new features ...done.\n",
      "[Sparsification] Loop: 0. Current number of features: 3\n",
      "TRAINING COMPLETE\n"
     ]
    }
   ],
   "source": [
    "# ---------------- INSTANTIATE WORKFLOW --------------------------------------\n",
    "training_data = generate_training_data()\n",
    "\n",
    "workflow = wf.WorkflowAF(feateng_steps = FEATENG_STEPS,\n",
    "                        units =  UNITS,\n",
    "                        featsel_runs = FEATSEL_RUNS,\n",
    "                        transformations = TRANSFORMATIONS,\n",
    "                        xtrain = training_data[['c','T']], \n",
    "                        ytrain = training_data['k'], \n",
    "                        scaling_type = SCALING_TYPE,\n",
    "                        stdalpha =  STD_ALPHA, \n",
    "                        rejection_thresshold = REJECTION_THR, \n",
    "                        fit_intercept = FIT_INTERCEPT) \n",
    "\n",
    "# ---------------- RUN TRAINING WORKFLOW --------------------------------------\n",
    "trained_workflow = workflow.run_workflow()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The trained model, its predictions and evaluation\n",
    "The result of the run_workflow() method is a TrainedWorkflow object, which stores the resulting coefficients after training, the expression found and has a handy method predict() to apply the model to new data. We use the predict method to compare the real target with the predction from the model, and compute regression metrics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean</th>\n",
       "      <th>stdev</th>\n",
       "      <th>coeff</th>\n",
       "      <th>coeff stdev</th>\n",
       "      <th>coeff |t|</th>\n",
       "      <th>coeff_corr</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>c</th>\n",
       "      <td>1.550000</td>\n",
       "      <td>0.854075</td>\n",
       "      <td>11.330666</td>\n",
       "      <td>0.001046</td>\n",
       "      <td>10834.582018</td>\n",
       "      <td>13.266590</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>T**(3/2)*sqrt(T*c)</th>\n",
       "      <td>97951.021798</td>\n",
       "      <td>36756.952537</td>\n",
       "      <td>2.445561</td>\n",
       "      <td>0.001046</td>\n",
       "      <td>2338.488253</td>\n",
       "      <td>0.000067</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>T</th>\n",
       "      <td>286.500000</td>\n",
       "      <td>25.622259</td>\n",
       "      <td>-0.186925</td>\n",
       "      <td>0.001046</td>\n",
       "      <td>178.741006</td>\n",
       "      <td>-0.007295</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                            mean         stdev      coeff  coeff stdev  \\\n",
       "c                       1.550000      0.854075  11.330666     0.001046   \n",
       "T**(3/2)*sqrt(T*c)  97951.021798  36756.952537   2.445561     0.001046   \n",
       "T                     286.500000     25.622259  -0.186925     0.001046   \n",
       "\n",
       "                       coeff |t|  coeff_corr  \n",
       "c                   10834.582018   13.266590  \n",
       "T**(3/2)*sqrt(T*c)   2338.488253    0.000067  \n",
       "T                     178.741006   -0.007295  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# ------------------------- INSPECT SYMBOLIC MODEL --------------------------\n",
    "trained_workflow.coeff_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 6.653328793234911 \\cdot 10^{-5} T^{\\frac{3}{2}} \\sqrt{T c} - 0.0072954162164353764 T + 13.266589527514514 c$"
      ],
      "text/plain": [
       "6.653328793234911e-5*T**(3/2)*sqrt(T*c) - 0.0072954162164353764*T + 13.266589527514514*c"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trained_workflow.eqn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mse: 0.1367177906453799\n",
      "r2: 0.9992615352216608\n",
      "mape: 0.02601023550725349\n"
     ]
    }
   ],
   "source": [
    "# ------------------------- EVALUATE SYMBOLIC MODEL --------------------------\n",
    "import sklearn.metrics as skmetrics \n",
    "\n",
    "y_real = training_data['k']\n",
    "y_hat = trained_workflow.predict(x = training_data[['c','T']])\n",
    "\n",
    "print('mse: {}'.format(skmetrics.mean_squared_error(y_real, y_hat)))\n",
    "print('r2: {}'.format(skmetrics.r2_score(y_real, y_hat)))\n",
    "print('mape: {}'.format(skmetrics.mean_absolute_percentage_error(y_real, y_hat)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "count    125000.000000\n",
       "mean         24.972218\n",
       "std          13.606592\n",
       "min           1.543477\n",
       "25%          13.230266\n",
       "50%          25.054764\n",
       "75%          36.672468\n",
       "max          51.068246\n",
       "Name: k, dtype: float64"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "training_data['k'].describe()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The R-squared coefficient indicates that almost all the variance in the (synthetic) target *k* can be explained by the symbolic model. Moreover, the resulting MSE is very low relative to the range and standard deviation of the target values.\n",
    "\n",
    "**Note:** The resulting equation is not the same that generated the data, but is simple (only 3 terms) and is a good surrogate to the initial equation given the low mse and high r2 metrics. As described in the manuscript, slight changes in data migh result in different expressions equally accurate. Model selection can be improved by, e.g. training in multiple subsamples of the training set and identifying consistent models, or implementing constraints."
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "1baa965d5efe3ac65b79dfc60c0d706280b1da80fedb7760faf2759126c4f253"
  },
  "kernelspec": {
   "display_name": "Python 3.8.6 64-bit",
   "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.1"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}