SymScore.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "676438a4",
   "metadata": {},
   "source": [
    "# SymScore (Symbolic Regression-Based Clinical Score Generator)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "dbc87e1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import re, copy, time\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import *\n",
    "from sklearn.metrics import confusion_matrix\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "from sklearn.cluster import KMeans\n",
    "\n",
    "import sklearn.metrics as metrics\n",
    "\n",
    "from tqdm.notebook import tqdm\n",
    "from gplearn.functions import make_function\n",
    "from gplearn.fitness import make_fitness\n",
    "from gplearn.genetic import SymbolicRegressor as SR\n",
    "from gplearn.genetic import SymbolicClassifier as SC\n",
    "from sympy import sympify,sin, cos\n",
    "from sklearn.metrics import mean_absolute_error\n",
    "from sklearn.model_selection import cross_val_score, KFold\n",
    "from sklearn.model_selection import StratifiedKFold, cross_val_predict\n",
    "from skopt import gp_minimize\n",
    "from skopt.space import Real, Integer\n",
    "import itertools\n",
    "import re\n",
    "import pickle"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b24b30b",
   "metadata": {},
   "source": [
    "# Input the excel file containing survey response\n",
    "\n",
    "Here, simply input the name of the Excel file containing the survey response. For example:\n",
    "\n",
    "filename = \"survery_response.xlsx\"\n",
    "\n",
    "The Excel file must contain the input variables as column names (e.g age, weight, BMI, etc) and the target variables. For regression tasks, the target variable is the total score of the input variables and for the classification tasks, the target variable is the diagnosis, where the entry is 0 for not having the disease and 1 as not having the disease. \n",
    "\n",
    "Furthermore, the specific task should be specified here. For the variable task:\n",
    "\n",
    "Write 0 for regression task and write 1 for classification task. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "b833ec26",
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = \"MCQI score only_kor_원본.xlsx\"\n",
    "task = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e55bc9b",
   "metadata": {},
   "source": [
    "# Inputs for regression task\n",
    "\n",
    "(Ignore this if your task is a classification task)\n",
    "Here, please input the necessary information regarding the shortened questionnaire responses.\n",
    "\n",
    "1. Selected categorical variables\n",
    "    \n",
    "    For categorical questions where the responses increase with increasing severity, please provide the column names for the selected categorical variables and the response list.\n",
    "    \n",
    "    For example, if the selected questions have column names 23, 28, 39, 51, 58, 60, then write selected_categorical_variables = [23, 28, 39, 51, 58, 60]\n",
    "    \n",
    "2. Categorical response list\n",
    "    \n",
    "    If the possible responses are 1 - No, 2 - Slight, 3 - Moderate, 4 - Severe, then write down\n",
    "categorical_response_list = [1, 2, 3, 4]\n",
    "\n",
    "3. Maximum number of items in questionnaire\n",
    "\n",
    "    Please also input the maximum number of items in the questionnaire. For example, for the MCQI, there are 60 total items so max_items = 60.\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "947a25c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Inputs for regression task\n",
    "if task == 0:\n",
    "    selected_categorical_variables = [23, 28, 39, 51, 58, 60]\n",
    "    categorical_response_list = [1, 2, 3, 4]\n",
    "    max_items = 60"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1db8743c",
   "metadata": {},
   "source": [
    "# Inputs for classification task\n",
    "\n",
    "(Ignore this if your task is a regression task) \n",
    "Here, please input the necessary information regarding the shortened questionnaire responses.\n",
    "\n",
    "1. Diagnosis column \n",
    "\n",
    "    Here, indicate the column name that contains the target value. In this case, the diagnosis for COMISA has column name 1. OSA so we write diagnosis_column = \"1. OSA\".\n",
    "    \n",
    "2. Selected categorical variables \n",
    "\n",
    "    For categorical questions where the responses increase with increasing severity, please provide the column names for the selected categorical variables and the response list.\n",
    "    For example, if the selected questions have column names \"ISI1a\", \"ISI1b\", \"ISI1c\", \"ISI2\", \"ISI5\", then write categorical_items = [\"ISI1a\", \"ISI1b\", \"ISI1c\", \"ISI2\", \"ISI5\"].\n",
    "    \n",
    "    \n",
    "    \n",
    "3. Categorical response list\n",
    "    \n",
    "    Additionally, if the possible responses are 1 - No, 2 - Slight, 3 - Moderate, 4 - Severe, then write down categorical_response_list = [1, 2, 3, 4]\n",
    "    \n",
    "4. Binary variables\n",
    "\n",
    "    If the questions contain binary variables, then indicate it. \n",
    "    For example, for the variable sex, write\n",
    "    categorical_items_binary = ['sex'] and \n",
    "    binary_variables = ['Male', 'Female']\n",
    "    \n",
    "5. Float variablles\n",
    "\n",
    "    Similarly indicate the float variables by writing the column names. \n",
    "    For example here, float_items = [\"age\", \"weight\", \"BMI\"].\n",
    "    \n",
    "6. Step size \n",
    "\n",
    "    Please also indicate the step sizes for the float items. These step sizes indicate how accurate you want your response grouping to be. For example, a step size of 1 will check if the grouping is optimal for every 1-step increase in value, while a step size of 10 will check for groupings after every 10 unit increase.\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "59922832",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Inputs for classification task\n",
    "if task == 1:\n",
    "    diagnosis_column = \"9. COMISA\"\n",
    "    categorical_items = [\"ISI1a\", \"ISI1b\", \"ISI1c\", \"ISI2\", \"ISI5\"]\n",
    "    categorical_response_list = [0, 1, 2, 3, 4]\n",
    "    categorical_items_binary = ['sex']\n",
    "    binary_variables = [\"Male\", \"Female\"]\n",
    "    float_items = [\"age\", \"weight\", \"BMI\"]\n",
    "    user_defined_step_sizes = {\n",
    "            'age': 10,\n",
    "            'weight': 10,\n",
    "            'BMI': 5\n",
    "        }\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea7b1d2f",
   "metadata": {},
   "source": [
    "# Start of SymScore framework"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "a23d7aad",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df = pd.read_excel(filename) \n",
    "df= df.replace(999,-1) #removing missing values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "30f088cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Necessary inputs for regression task\n",
    "if task == 0:\n",
    "    categorical_items = [k for k in range(1,max_items+1)] \n",
    "    df_pre = df[selected_categorical_variables]\n",
    "    df_linear = df[selected_categorical_variables]\n",
    "    max_total_response = max_items*max(categorical_response_list)\n",
    "\n",
    "elif task ==1:\n",
    "    max_index = len(df[df[diagnosis_column] == 1])\n",
    "\n",
    "    #Combine all items\n",
    "    all_items = categorical_items + categorical_items_binary + float_items\n",
    "\n",
    "    # Determine the maximum value for each float item\n",
    "    max_values = {item: max(df[item]) for item in float_items}\n",
    "\n",
    "    # Generate the range parameters dynamically based on the max values and user-defined step sizes\n",
    "    range_params = {item: (0, max_values[item], user_defined_step_sizes[item]) for item in float_items}\n",
    "\n",
    "    # Generate the cutoffs dictionary dynamically\n",
    "    cutoffs = {item: np.arange(*range_params[item]) for item in float_items}\n",
    "    total_variable_length = len(categorical_items+float_items)\n",
    "    df_neg = df.iloc[df[df[diagnosis_column] == 0].index[:len(df[df[diagnosis_column] == 1])]]\n",
    "    df_pos = df.iloc[df[df[diagnosis_column] == 1].index[:len(df[df[diagnosis_column] == 1])]]\n",
    "    df_linear = pd.concat([df_pos,df_neg], axis=0)\n",
    "    df_pre = df_linear[categorical_items+float_items+categorical_items_binary]    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "b0559585",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/zn/b4jwzml16gq9cw9j3gdvbwt80000gn/T/ipykernel_10634/2824330735.py:8: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  df_pre[\"Thr{}(I{})\".format(j,k)] = df[k][1:].apply(lambda x: 1 if x>=j else 0)\n",
      "/var/folders/zn/b4jwzml16gq9cw9j3gdvbwt80000gn/T/ipykernel_10634/2824330735.py:8: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  df_pre[\"Thr{}(I{})\".format(j,k)] = df[k][1:].apply(lambda x: 1 if x>=j else 0)\n",
      "/var/folders/zn/b4jwzml16gq9cw9j3gdvbwt80000gn/T/ipykernel_10634/2824330735.py:8: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  df_pre[\"Thr{}(I{})\".format(j,k)] = df[k][1:].apply(lambda x: 1 if x>=j else 0)\n",
      "/var/folders/zn/b4jwzml16gq9cw9j3gdvbwt80000gn/T/ipykernel_10634/2824330735.py:8: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  df_pre[\"Thr{}(I{})\".format(j,k)] = df[k][1:].apply(lambda x: 1 if x>=j else 0)\n",
      "/var/folders/zn/b4jwzml16gq9cw9j3gdvbwt80000gn/T/ipykernel_10634/2824330735.py:8: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  df_pre[\"Thr{}(I{})\".format(j,k)] = df[k][1:].apply(lambda x: 1 if x>=j else 0)\n",
      "/var/folders/zn/b4jwzml16gq9cw9j3gdvbwt80000gn/T/ipykernel_10634/2824330735.py:8: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  df_pre[\"Thr{}(I{})\".format(j,k)] = df[k][1:].apply(lambda x: 1 if x>=j else 0)\n",
      "/var/folders/zn/b4jwzml16gq9cw9j3gdvbwt80000gn/T/ipykernel_10634/2824330735.py:8: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  df_pre[\"Thr{}(I{})\".format(j,k)] = df[k][1:].apply(lambda x: 1 if x>=j else 0)\n",
      "/var/folders/zn/b4jwzml16gq9cw9j3gdvbwt80000gn/T/ipykernel_10634/2824330735.py:8: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  df_pre[\"Thr{}(I{})\".format(j,k)] = df[k][1:].apply(lambda x: 1 if x>=j else 0)\n",
      "/var/folders/zn/b4jwzml16gq9cw9j3gdvbwt80000gn/T/ipykernel_10634/2824330735.py:8: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  df_pre[\"Thr{}(I{})\".format(j,k)] = df[k][1:].apply(lambda x: 1 if x>=j else 0)\n",
      "/var/folders/zn/b4jwzml16gq9cw9j3gdvbwt80000gn/T/ipykernel_10634/2824330735.py:8: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  df_pre[\"Thr{}(I{})\".format(j,k)] = df[k][1:].apply(lambda x: 1 if x>=j else 0)\n",
      "/var/folders/zn/b4jwzml16gq9cw9j3gdvbwt80000gn/T/ipykernel_10634/2824330735.py:8: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  df_pre[\"Thr{}(I{})\".format(j,k)] = df[k][1:].apply(lambda x: 1 if x>=j else 0)\n",
      "/var/folders/zn/b4jwzml16gq9cw9j3gdvbwt80000gn/T/ipykernel_10634/2824330735.py:8: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  df_pre[\"Thr{}(I{})\".format(j,k)] = df[k][1:].apply(lambda x: 1 if x>=j else 0)\n",
      "/var/folders/zn/b4jwzml16gq9cw9j3gdvbwt80000gn/T/ipykernel_10634/2824330735.py:8: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  df_pre[\"Thr{}(I{})\".format(j,k)] = df[k][1:].apply(lambda x: 1 if x>=j else 0)\n",
      "/var/folders/zn/b4jwzml16gq9cw9j3gdvbwt80000gn/T/ipykernel_10634/2824330735.py:8: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  df_pre[\"Thr{}(I{})\".format(j,k)] = df[k][1:].apply(lambda x: 1 if x>=j else 0)\n",
      "/var/folders/zn/b4jwzml16gq9cw9j3gdvbwt80000gn/T/ipykernel_10634/2824330735.py:8: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  df_pre[\"Thr{}(I{})\".format(j,k)] = df[k][1:].apply(lambda x: 1 if x>=j else 0)\n",
      "/var/folders/zn/b4jwzml16gq9cw9j3gdvbwt80000gn/T/ipykernel_10634/2824330735.py:8: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  df_pre[\"Thr{}(I{})\".format(j,k)] = df[k][1:].apply(lambda x: 1 if x>=j else 0)\n",
      "/var/folders/zn/b4jwzml16gq9cw9j3gdvbwt80000gn/T/ipykernel_10634/2824330735.py:8: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  df_pre[\"Thr{}(I{})\".format(j,k)] = df[k][1:].apply(lambda x: 1 if x>=j else 0)\n",
      "/var/folders/zn/b4jwzml16gq9cw9j3gdvbwt80000gn/T/ipykernel_10634/2824330735.py:8: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  df_pre[\"Thr{}(I{})\".format(j,k)] = df[k][1:].apply(lambda x: 1 if x>=j else 0)\n",
      "/var/folders/zn/b4jwzml16gq9cw9j3gdvbwt80000gn/T/ipykernel_10634/2824330735.py:9: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  df_pre.drop(columns = selected_categorical_variables, inplace = True)\n"
     ]
    }
   ],
   "source": [
    "# Creates new binary indicator columns in the DataFrame df_pre, \n",
    "# based on whether values in specific columns of the DataFrame df meet certain threshold conditions. \n",
    "# The thresholds are defined by the range of integers in the inner loop (j).\n",
    "\n",
    "if task == 0:\n",
    "    for k in selected_categorical_variables:\n",
    "        for j in range(min(categorical_response_list), len(categorical_response_list)):\n",
    "            df_pre[\"Thr{}(I{})\".format(j,k)] = df[k][1:].apply(lambda x: 1 if x>=j else 0)\n",
    "    df_pre.drop(columns = selected_categorical_variables, inplace = True)\n",
    "    \n",
    "elif task ==1:\n",
    "    for k in categorical_items:\n",
    "        for j in range(min(categorical_response_list), len(categorical_response_list)+1):\n",
    "            df_pre[\"Thr{}({})\".format(j,k)] = df_pre[k].apply(lambda x: 1 if x>=j else 0)\n",
    "            df_pre[\"iThr{}({})\".format(j,k)] = df_pre[k].apply(lambda x: 1 if x<j else 0)\n",
    "\n",
    "    df_pre.drop(columns = categorical_items, inplace = True)\n",
    "\n",
    "    \n",
    "    # Generalized loop to create threshold columns\n",
    "    for var, cut_values in cutoffs.items():\n",
    "        for j in cut_values:\n",
    "            # Formatting the column names with properly placed decimal points\n",
    "            col_name_thr = f\"Thr{var.capitalize()}{j:.0f}({var})\"\n",
    "            col_name_i_thr = f\"iThr{var.capitalize()}{j:.0f}({var})\"\n",
    "\n",
    "            # Create threshold columns\n",
    "            df_pre[col_name_thr] = df_pre[var].apply(lambda x: 1 if x >= j else 0)\n",
    "            df_pre[col_name_i_thr] = df_pre[var].apply(lambda x: 1 if x < j else 0)\n",
    "\n",
    "    \n",
    "    df_pre[binary_variables[0]] = df_pre[categorical_items_binary[0]].apply(lambda x: 1 if x==0 else 0)\n",
    "    df_pre[binary_variables[1]] = df_pre[categorical_items_binary[0]].apply(lambda x: 1 if x==1 else 0)\n",
    "\n",
    "\n",
    "    df_pre.drop(columns = float_items, inplace = True)\n",
    "    df_pre.drop(columns = categorical_items_binary, inplace = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99301a31",
   "metadata": {},
   "source": [
    "# Prepare the data for training \n",
    "Here the data is split into training and testing sets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "1c191337",
   "metadata": {},
   "outputs": [],
   "source": [
    "if task ==0:\n",
    "    X = df_pre[1:] # predictor variable\n",
    "    y = df[categorical_items][:][1:].sum(axis=1).values # target variable\n",
    "    x_train, x_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=1) # split the predictor and target variables (X and y) into training and testing sets \n",
    "\n",
    "    alpha = 13 # number of synthetic data points to be added\n",
    "    added_data_a = dict([(keyname, [1]*alpha) for keyname in x_train.columns])\n",
    "    tablea = pd.DataFrame(added_data_a)\n",
    "    x_train = pd.concat([x_train, tablea])\n",
    "\n",
    "    for i in range(alpha):\n",
    "        y_train = np.append(y_train, [max_total_response])\n",
    "    \n",
    "        \n",
    "elif task==1:\n",
    "    X = df_pre/total_variable_length \n",
    "    y = df_linear[diagnosis_column]\n",
    "    x_train, x_test, y_train, y_test = train_test_split(X, y, stratify=y, test_size=0.3, random_state=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c728a6a",
   "metadata": {},
   "source": [
    "# Training \n",
    "\n",
    "Initialize and train a symbolic regression model using only the addition function, which means the model will attempt to construct mathematical expressions using only addition operations to fit the training data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "4a4c2bdc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the objective function for regression task\n",
    "def sym_regression(params):\n",
    "    p_subtree_mutation, p_hoist_mutation, p_point_mutation = params\n",
    "    p_crossover=0.7\n",
    "    # Penalty\n",
    "    # Ensure the sum of probabilities is <= 1.0\n",
    "    if p_crossover + p_subtree_mutation + p_hoist_mutation + p_point_mutation > 1.0:\n",
    "        return 1e6  \n",
    "    \n",
    "    est_gp = SR(\n",
    "        population_size=10000,\n",
    "        function_set=function_set,\n",
    "        generations=20,\n",
    "        stopping_criteria=0.01,\n",
    "        p_crossover=0.7,\n",
    "        p_subtree_mutation=p_subtree_mutation,\n",
    "        p_hoist_mutation=p_hoist_mutation,\n",
    "        p_point_mutation=p_point_mutation,\n",
    "        max_samples=0.9,\n",
    "        verbose=0,\n",
    "        metric=\"mean absolute error\",\n",
    "        parsimony_coefficient=0.01,\n",
    "        random_state=10,\n",
    "        feature_names=x_train.columns\n",
    "    )\n",
    "    \n",
    "    est_gp.fit(x_train, y_train)\n",
    "    return -est_gp.score(x_test, y_test)\n",
    "\n",
    "# Define the objective function for classification task\n",
    "def sym_classification(params):\n",
    "    p_subtree_mutation, p_hoist_mutation, p_point_mutation = params\n",
    "    p_crossover=0.7\n",
    "    # Penalty\n",
    "    # Ensure the sum of probabilities is <= 1.0\n",
    "    if p_crossover + p_subtree_mutation + p_hoist_mutation + p_point_mutation > 1.0:\n",
    "        return 1e6\n",
    "    \n",
    "    est_gp = SC(\n",
    "        population_size=20000,\n",
    "        function_set=function_set,\n",
    "        generations=50,\n",
    "        stopping_criteria=10**(-4),\n",
    "        p_crossover=0.7,\n",
    "        p_subtree_mutation=p_subtree_mutation,\n",
    "        p_hoist_mutation=p_hoist_mutation,\n",
    "        p_point_mutation=p_point_mutation,\n",
    "        max_samples=0.9,\n",
    "        verbose=0,\n",
    "        random_state=1,\n",
    "        feature_names=x_train.columns\n",
    "    )\n",
    "    \n",
    "    est_gp.fit(x_train, y_train)\n",
    "    predictions_proba = est_gp.predict_proba(x_test)[:, 1]  # Get probabilities for the positive class\n",
    "    return -roc_auc_score(y_test, predictions_proba)\n",
    "\n",
    "# Define the search space for Task 0\n",
    "pbounds_regression = [\n",
    "    Real(0.0, 0.1, name='p_hoist_mutation'),\n",
    "    Real(0.0, 0.1, name='p_subtree_mutation'),\n",
    "    Real(0.0, 0.1, name='p_point_mutation'),\n",
    "]\n",
    "\n",
    "# Define the search space for Task 1\n",
    "pbounds_classification = [\n",
    "    Real(0.0, 0.1, name='p_hoist_mutation'),\n",
    "    Real(0.0, 0.1, name='p_subtree_mutation'),\n",
    "    Real(0.0, 0.1, name='p_point_mutation'),\n",
    "]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "b38d187d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration No: 1 started. Evaluating function at random point.\n",
      "Iteration No: 1 ended. Evaluation done at random point.\n",
      "Time taken: 140.6171\n",
      "Function value obtained: -0.8190\n",
      "Current minimum: -0.8190\n",
      "Iteration No: 2 started. Evaluating function at random point.\n",
      "Iteration No: 2 ended. Evaluation done at random point.\n",
      "Time taken: 143.0240\n",
      "Function value obtained: -0.8174\n",
      "Current minimum: -0.8190\n",
      "Iteration No: 3 started. Evaluating function at random point.\n",
      "Iteration No: 3 ended. Evaluation done at random point.\n",
      "Time taken: 142.0354\n",
      "Function value obtained: -0.8044\n",
      "Current minimum: -0.8190\n",
      "Iteration No: 4 started. Evaluating function at random point.\n",
      "Iteration No: 4 ended. Evaluation done at random point.\n",
      "Time taken: 140.2379\n",
      "Function value obtained: -0.8203\n",
      "Current minimum: -0.8203\n",
      "Iteration No: 5 started. Evaluating function at random point.\n",
      "Iteration No: 5 ended. Evaluation done at random point.\n",
      "Time taken: 138.7180\n",
      "Function value obtained: -0.8133\n",
      "Current minimum: -0.8203\n",
      "Iteration No: 6 started. Evaluating function at random point.\n",
      "Iteration No: 6 ended. Evaluation done at random point.\n",
      "Time taken: 142.0296\n",
      "Function value obtained: -0.8079\n",
      "Current minimum: -0.8203\n",
      "Iteration No: 7 started. Evaluating function at random point.\n",
      "Iteration No: 7 ended. Evaluation done at random point.\n",
      "Time taken: 140.2385\n",
      "Function value obtained: -0.8164\n",
      "Current minimum: -0.8203\n",
      "Iteration No: 8 started. Evaluating function at random point.\n",
      "Iteration No: 8 ended. Evaluation done at random point.\n",
      "Time taken: 142.1763\n",
      "Function value obtained: -0.8096\n",
      "Current minimum: -0.8203\n",
      "Iteration No: 9 started. Evaluating function at random point.\n",
      "Iteration No: 9 ended. Evaluation done at random point.\n",
      "Time taken: 140.4530\n",
      "Function value obtained: -0.8120\n",
      "Current minimum: -0.8203\n",
      "Iteration No: 10 started. Evaluating function at random point.\n",
      "Iteration No: 10 ended. Evaluation done at random point.\n",
      "Time taken: 139.8207\n",
      "Function value obtained: -0.8196\n",
      "Current minimum: -0.8203\n",
      "Iteration No: 11 started. Searching for the next optimal point.\n",
      "Iteration No: 11 ended. Search finished for the next optimal point.\n",
      "Time taken: 144.9131\n",
      "Function value obtained: -0.8098\n",
      "Current minimum: -0.8203\n",
      "Iteration No: 12 started. Searching for the next optimal point.\n",
      "Iteration No: 12 ended. Search finished for the next optimal point.\n",
      "Time taken: 145.5584\n",
      "Function value obtained: -0.8096\n",
      "Current minimum: -0.8203\n",
      "Iteration No: 13 started. Searching for the next optimal point.\n",
      "Iteration No: 13 ended. Search finished for the next optimal point.\n",
      "Time taken: 143.3649\n",
      "Function value obtained: -0.8203\n",
      "Current minimum: -0.8203\n",
      "Iteration No: 14 started. Searching for the next optimal point.\n",
      "Iteration No: 14 ended. Search finished for the next optimal point.\n",
      "Time taken: 140.9114\n",
      "Function value obtained: -0.8122\n",
      "Current minimum: -0.8203\n",
      "Iteration No: 15 started. Searching for the next optimal point.\n",
      "Iteration No: 15 ended. Search finished for the next optimal point.\n",
      "Time taken: 138.3079\n",
      "Function value obtained: -0.8187\n",
      "Current minimum: -0.8203\n",
      "Iteration No: 16 started. Searching for the next optimal point.\n",
      "Iteration No: 16 ended. Search finished for the next optimal point.\n",
      "Time taken: 136.2421\n",
      "Function value obtained: -0.8196\n",
      "Current minimum: -0.8203\n",
      "Iteration No: 17 started. Searching for the next optimal point.\n",
      "Iteration No: 17 ended. Search finished for the next optimal point.\n",
      "Time taken: 136.7301\n",
      "Function value obtained: -0.8163\n",
      "Current minimum: -0.8203\n",
      "Iteration No: 18 started. Searching for the next optimal point.\n",
      "Iteration No: 18 ended. Search finished for the next optimal point.\n",
      "Time taken: 139.5965\n",
      "Function value obtained: -0.8136\n",
      "Current minimum: -0.8203\n",
      "Iteration No: 19 started. Searching for the next optimal point.\n",
      "Iteration No: 19 ended. Search finished for the next optimal point.\n",
      "Time taken: 137.9848\n",
      "Function value obtained: -0.8184\n",
      "Current minimum: -0.8203\n",
      "Iteration No: 20 started. Searching for the next optimal point.\n",
      "Iteration No: 20 ended. Search finished for the next optimal point.\n",
      "Time taken: 138.7904\n",
      "Function value obtained: -0.8163\n",
      "Current minimum: -0.8203\n",
      "Best parameters for Task 0: [0.08463109166860174, 0.03132735169322752, 0.052454815957287154]\n",
      "    |   Population Average    |             Best Individual              |\n",
      "---- ------------------------- ------------------------------------------ ----------\n",
      " Gen   Length          Fitness   Length          Fitness      OOB Fitness  Time Left\n",
      "   0    26.19          132.725      127          90.4492          84.8357     53.46s\n",
      "   1   118.61           100.31      253          47.9677          44.2044      1.22m\n",
      "   2   166.06           82.871      379          14.0698          18.7604      1.60m\n",
      "   3   264.68          49.8119      379           13.909           15.637      1.86m\n",
      "   4   366.25          23.9586      367          13.4216           18.188      1.97m\n",
      "   5   377.21          19.7194      387          12.9329           22.961      1.81m\n",
      "   6   372.45          20.1126      383          12.9229          23.7089      1.99m\n",
      "   7   368.07           19.952      363          12.6451           24.804      1.57m\n",
      "   8   366.12          20.0124      359          12.5635          27.6277      1.47m\n",
      "   9   364.11          20.1106      365          12.6972          23.9194      1.27m\n",
      "  10   364.23          19.9407      371          12.2806          26.1752      1.24m\n",
      "  11   363.66          20.0496      357          12.5994           25.562     59.61s\n",
      "  12   364.45          20.1884      369          12.3293          28.2967     52.52s\n",
      "  13   365.16          20.0392      373          12.4432          22.6701     44.81s\n",
      "  14   364.95          19.4797      369           12.342          24.1548     37.46s\n",
      "  15   365.49          19.4112      369          12.3764          24.2006     31.77s\n",
      "  16   364.92          19.2989      367           12.044          26.8224     22.45s\n",
      "  17   363.98          19.4111      363          12.0733          27.8905     14.99s\n",
      "  18   365.57          19.1954      365          12.2773          24.3032      7.50s\n",
      "  19   364.50          19.1075      371          12.1918          26.3289      0.00s\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<style>#sk-container-id-2 {\n",
       "  /* Definition of color scheme common for light and dark mode */\n",
       "  --sklearn-color-text: black;\n",
       "  --sklearn-color-line: gray;\n",
       "  /* Definition of color scheme for unfitted estimators */\n",
       "  --sklearn-color-unfitted-level-0: #fff5e6;\n",
       "  --sklearn-color-unfitted-level-1: #f6e4d2;\n",
       "  --sklearn-color-unfitted-level-2: #ffe0b3;\n",
       "  --sklearn-color-unfitted-level-3: chocolate;\n",
       "  /* Definition of color scheme for fitted estimators */\n",
       "  --sklearn-color-fitted-level-0: #f0f8ff;\n",
       "  --sklearn-color-fitted-level-1: #d4ebff;\n",
       "  --sklearn-color-fitted-level-2: #b3dbfd;\n",
       "  --sklearn-color-fitted-level-3: cornflowerblue;\n",
       "\n",
       "  /* Specific color for light theme */\n",
       "  --sklearn-color-text-on-default-background: var(--sg-text-color, var(--theme-code-foreground, var(--jp-content-font-color1, black)));\n",
       "  --sklearn-color-background: var(--sg-background-color, var(--theme-background, var(--jp-layout-color0, white)));\n",
       "  --sklearn-color-border-box: var(--sg-text-color, var(--theme-code-foreground, var(--jp-content-font-color1, black)));\n",
       "  --sklearn-color-icon: #696969;\n",
       "\n",
       "  @media (prefers-color-scheme: dark) {\n",
       "    /* Redefinition of color scheme for dark theme */\n",
       "    --sklearn-color-text-on-default-background: var(--sg-text-color, var(--theme-code-foreground, var(--jp-content-font-color1, white)));\n",
       "    --sklearn-color-background: var(--sg-background-color, var(--theme-background, var(--jp-layout-color0, #111)));\n",
       "    --sklearn-color-border-box: var(--sg-text-color, var(--theme-code-foreground, var(--jp-content-font-color1, white)));\n",
       "    --sklearn-color-icon: #878787;\n",
       "  }\n",
       "}\n",
       "\n",
       "#sk-container-id-2 {\n",
       "  color: var(--sklearn-color-text);\n",
       "}\n",
       "\n",
       "#sk-container-id-2 pre {\n",
       "  padding: 0;\n",
       "}\n",
       "\n",
       "#sk-container-id-2 input.sk-hidden--visually {\n",
       "  border: 0;\n",
       "  clip: rect(1px 1px 1px 1px);\n",
       "  clip: rect(1px, 1px, 1px, 1px);\n",
       "  height: 1px;\n",
       "  margin: -1px;\n",
       "  overflow: hidden;\n",
       "  padding: 0;\n",
       "  position: absolute;\n",
       "  width: 1px;\n",
       "}\n",
       "\n",
       "#sk-container-id-2 div.sk-dashed-wrapped {\n",
       "  border: 1px dashed var(--sklearn-color-line);\n",
       "  margin: 0 0.4em 0.5em 0.4em;\n",
       "  box-sizing: border-box;\n",
       "  padding-bottom: 0.4em;\n",
       "  background-color: var(--sklearn-color-background);\n",
       "}\n",
       "\n",
       "#sk-container-id-2 div.sk-container {\n",
       "  /* jupyter's `normalize.less` sets `[hidden] { display: none; }`\n",
       "     but bootstrap.min.css set `[hidden] { display: none !important; }`\n",
       "     so we also need the `!important` here to be able to override the\n",
       "     default hidden behavior on the sphinx rendered scikit-learn.org.\n",
       "     See: /~https://github.com/scikit-learn/scikit-learn/issues/21755 */\n",
       "  display: inline-block !important;\n",
       "  position: relative;\n",
       "}\n",
       "\n",
       "#sk-container-id-2 div.sk-text-repr-fallback {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       "div.sk-parallel-item,\n",
       "div.sk-serial,\n",
       "div.sk-item {\n",
       "  /* draw centered vertical line to link estimators */\n",
       "  background-image: linear-gradient(var(--sklearn-color-text-on-default-background), var(--sklearn-color-text-on-default-background));\n",
       "  background-size: 2px 100%;\n",
       "  background-repeat: no-repeat;\n",
       "  background-position: center center;\n",
       "}\n",
       "\n",
       "/* Parallel-specific style estimator block */\n",
       "\n",
       "#sk-container-id-2 div.sk-parallel-item::after {\n",
       "  content: \"\";\n",
       "  width: 100%;\n",
       "  border-bottom: 2px solid var(--sklearn-color-text-on-default-background);\n",
       "  flex-grow: 1;\n",
       "}\n",
       "\n",
       "#sk-container-id-2 div.sk-parallel {\n",
       "  display: flex;\n",
       "  align-items: stretch;\n",
       "  justify-content: center;\n",
       "  background-color: var(--sklearn-color-background);\n",
       "  position: relative;\n",
       "}\n",
       "\n",
       "#sk-container-id-2 div.sk-parallel-item {\n",
       "  display: flex;\n",
       "  flex-direction: column;\n",
       "}\n",
       "\n",
       "#sk-container-id-2 div.sk-parallel-item:first-child::after {\n",
       "  align-self: flex-end;\n",
       "  width: 50%;\n",
       "}\n",
       "\n",
       "#sk-container-id-2 div.sk-parallel-item:last-child::after {\n",
       "  align-self: flex-start;\n",
       "  width: 50%;\n",
       "}\n",
       "\n",
       "#sk-container-id-2 div.sk-parallel-item:only-child::after {\n",
       "  width: 0;\n",
       "}\n",
       "\n",
       "/* Serial-specific style estimator block */\n",
       "\n",
       "#sk-container-id-2 div.sk-serial {\n",
       "  display: flex;\n",
       "  flex-direction: column;\n",
       "  align-items: center;\n",
       "  background-color: var(--sklearn-color-background);\n",
       "  padding-right: 1em;\n",
       "  padding-left: 1em;\n",
       "}\n",
       "\n",
       "\n",
       "/* Toggleable style: style used for estimator/Pipeline/ColumnTransformer box that is\n",
       "clickable and can be expanded/collapsed.\n",
       "- Pipeline and ColumnTransformer use this feature and define the default style\n",
       "- Estimators will overwrite some part of the style using the `sk-estimator` class\n",
       "*/\n",
       "\n",
       "/* Pipeline and ColumnTransformer style (default) */\n",
       "\n",
       "#sk-container-id-2 div.sk-toggleable {\n",
       "  /* Default theme specific background. It is overwritten whether we have a\n",
       "  specific estimator or a Pipeline/ColumnTransformer */\n",
       "  background-color: var(--sklearn-color-background);\n",
       "}\n",
       "\n",
       "/* Toggleable label */\n",
       "#sk-container-id-2 label.sk-toggleable__label {\n",
       "  cursor: pointer;\n",
       "  display: block;\n",
       "  width: 100%;\n",
       "  margin-bottom: 0;\n",
       "  padding: 0.5em;\n",
       "  box-sizing: border-box;\n",
       "  text-align: center;\n",
       "}\n",
       "\n",
       "#sk-container-id-2 label.sk-toggleable__label-arrow:before {\n",
       "  /* Arrow on the left of the label */\n",
       "  content: \"▸\";\n",
       "  float: left;\n",
       "  margin-right: 0.25em;\n",
       "  color: var(--sklearn-color-icon);\n",
       "}\n",
       "\n",
       "#sk-container-id-2 label.sk-toggleable__label-arrow:hover:before {\n",
       "  color: var(--sklearn-color-text);\n",
       "}\n",
       "\n",
       "/* Toggleable content - dropdown */\n",
       "\n",
       "#sk-container-id-2 div.sk-toggleable__content {\n",
       "  max-height: 0;\n",
       "  max-width: 0;\n",
       "  overflow: hidden;\n",
       "  text-align: left;\n",
       "  /* unfitted */\n",
       "  background-color: var(--sklearn-color-unfitted-level-0);\n",
       "}\n",
       "\n",
       "#sk-container-id-2 div.sk-toggleable__content.fitted {\n",
       "  /* fitted */\n",
       "  background-color: var(--sklearn-color-fitted-level-0);\n",
       "}\n",
       "\n",
       "#sk-container-id-2 div.sk-toggleable__content pre {\n",
       "  margin: 0.2em;\n",
       "  border-radius: 0.25em;\n",
       "  color: var(--sklearn-color-text);\n",
       "  /* unfitted */\n",
       "  background-color: var(--sklearn-color-unfitted-level-0);\n",
       "}\n",
       "\n",
       "#sk-container-id-2 div.sk-toggleable__content.fitted pre {\n",
       "  /* unfitted */\n",
       "  background-color: var(--sklearn-color-fitted-level-0);\n",
       "}\n",
       "\n",
       "#sk-container-id-2 input.sk-toggleable__control:checked~div.sk-toggleable__content {\n",
       "  /* Expand drop-down */\n",
       "  max-height: 200px;\n",
       "  max-width: 100%;\n",
       "  overflow: auto;\n",
       "}\n",
       "\n",
       "#sk-container-id-2 input.sk-toggleable__control:checked~label.sk-toggleable__label-arrow:before {\n",
       "  content: \"▾\";\n",
       "}\n",
       "\n",
       "/* Pipeline/ColumnTransformer-specific style */\n",
       "\n",
       "#sk-container-id-2 div.sk-label input.sk-toggleable__control:checked~label.sk-toggleable__label {\n",
       "  color: var(--sklearn-color-text);\n",
       "  background-color: var(--sklearn-color-unfitted-level-2);\n",
       "}\n",
       "\n",
       "#sk-container-id-2 div.sk-label.fitted input.sk-toggleable__control:checked~label.sk-toggleable__label {\n",
       "  background-color: var(--sklearn-color-fitted-level-2);\n",
       "}\n",
       "\n",
       "/* Estimator-specific style */\n",
       "\n",
       "/* Colorize estimator box */\n",
       "#sk-container-id-2 div.sk-estimator input.sk-toggleable__control:checked~label.sk-toggleable__label {\n",
       "  /* unfitted */\n",
       "  background-color: var(--sklearn-color-unfitted-level-2);\n",
       "}\n",
       "\n",
       "#sk-container-id-2 div.sk-estimator.fitted input.sk-toggleable__control:checked~label.sk-toggleable__label {\n",
       "  /* fitted */\n",
       "  background-color: var(--sklearn-color-fitted-level-2);\n",
       "}\n",
       "\n",
       "#sk-container-id-2 div.sk-label label.sk-toggleable__label,\n",
       "#sk-container-id-2 div.sk-label label {\n",
       "  /* The background is the default theme color */\n",
       "  color: var(--sklearn-color-text-on-default-background);\n",
       "}\n",
       "\n",
       "/* On hover, darken the color of the background */\n",
       "#sk-container-id-2 div.sk-label:hover label.sk-toggleable__label {\n",
       "  color: var(--sklearn-color-text);\n",
       "  background-color: var(--sklearn-color-unfitted-level-2);\n",
       "}\n",
       "\n",
       "/* Label box, darken color on hover, fitted */\n",
       "#sk-container-id-2 div.sk-label.fitted:hover label.sk-toggleable__label.fitted {\n",
       "  color: var(--sklearn-color-text);\n",
       "  background-color: var(--sklearn-color-fitted-level-2);\n",
       "}\n",
       "\n",
       "/* Estimator label */\n",
       "\n",
       "#sk-container-id-2 div.sk-label label {\n",
       "  font-family: monospace;\n",
       "  font-weight: bold;\n",
       "  display: inline-block;\n",
       "  line-height: 1.2em;\n",
       "}\n",
       "\n",
       "#sk-container-id-2 div.sk-label-container {\n",
       "  text-align: center;\n",
       "}\n",
       "\n",
       "/* Estimator-specific */\n",
       "#sk-container-id-2 div.sk-estimator {\n",
       "  font-family: monospace;\n",
       "  border: 1px dotted var(--sklearn-color-border-box);\n",
       "  border-radius: 0.25em;\n",
       "  box-sizing: border-box;\n",
       "  margin-bottom: 0.5em;\n",
       "  /* unfitted */\n",
       "  background-color: var(--sklearn-color-unfitted-level-0);\n",
       "}\n",
       "\n",
       "#sk-container-id-2 div.sk-estimator.fitted {\n",
       "  /* fitted */\n",
       "  background-color: var(--sklearn-color-fitted-level-0);\n",
       "}\n",
       "\n",
       "/* on hover */\n",
       "#sk-container-id-2 div.sk-estimator:hover {\n",
       "  /* unfitted */\n",
       "  background-color: var(--sklearn-color-unfitted-level-2);\n",
       "}\n",
       "\n",
       "#sk-container-id-2 div.sk-estimator.fitted:hover {\n",
       "  /* fitted */\n",
       "  background-color: var(--sklearn-color-fitted-level-2);\n",
       "}\n",
       "\n",
       "/* Specification for estimator info (e.g. \"i\" and \"?\") */\n",
       "\n",
       "/* Common style for \"i\" and \"?\" */\n",
       "\n",
       ".sk-estimator-doc-link,\n",
       "a:link.sk-estimator-doc-link,\n",
       "a:visited.sk-estimator-doc-link {\n",
       "  float: right;\n",
       "  font-size: smaller;\n",
       "  line-height: 1em;\n",
       "  font-family: monospace;\n",
       "  background-color: var(--sklearn-color-background);\n",
       "  border-radius: 1em;\n",
       "  height: 1em;\n",
       "  width: 1em;\n",
       "  text-decoration: none !important;\n",
       "  margin-left: 1ex;\n",
       "  /* unfitted */\n",
       "  border: var(--sklearn-color-unfitted-level-1) 1pt solid;\n",
       "  color: var(--sklearn-color-unfitted-level-1);\n",
       "}\n",
       "\n",
       ".sk-estimator-doc-link.fitted,\n",
       "a:link.sk-estimator-doc-link.fitted,\n",
       "a:visited.sk-estimator-doc-link.fitted {\n",
       "  /* fitted */\n",
       "  border: var(--sklearn-color-fitted-level-1) 1pt solid;\n",
       "  color: var(--sklearn-color-fitted-level-1);\n",
       "}\n",
       "\n",
       "/* On hover */\n",
       "div.sk-estimator:hover .sk-estimator-doc-link:hover,\n",
       ".sk-estimator-doc-link:hover,\n",
       "div.sk-label-container:hover .sk-estimator-doc-link:hover,\n",
       ".sk-estimator-doc-link:hover {\n",
       "  /* unfitted */\n",
       "  background-color: var(--sklearn-color-unfitted-level-3);\n",
       "  color: var(--sklearn-color-background);\n",
       "  text-decoration: none;\n",
       "}\n",
       "\n",
       "div.sk-estimator.fitted:hover .sk-estimator-doc-link.fitted:hover,\n",
       ".sk-estimator-doc-link.fitted:hover,\n",
       "div.sk-label-container:hover .sk-estimator-doc-link.fitted:hover,\n",
       ".sk-estimator-doc-link.fitted:hover {\n",
       "  /* fitted */\n",
       "  background-color: var(--sklearn-color-fitted-level-3);\n",
       "  color: var(--sklearn-color-background);\n",
       "  text-decoration: none;\n",
       "}\n",
       "\n",
       "/* Span, style for the box shown on hovering the info icon */\n",
       ".sk-estimator-doc-link span {\n",
       "  display: none;\n",
       "  z-index: 9999;\n",
       "  position: relative;\n",
       "  font-weight: normal;\n",
       "  right: .2ex;\n",
       "  padding: .5ex;\n",
       "  margin: .5ex;\n",
       "  width: min-content;\n",
       "  min-width: 20ex;\n",
       "  max-width: 50ex;\n",
       "  color: var(--sklearn-color-text);\n",
       "  box-shadow: 2pt 2pt 4pt #999;\n",
       "  /* unfitted */\n",
       "  background: var(--sklearn-color-unfitted-level-0);\n",
       "  border: .5pt solid var(--sklearn-color-unfitted-level-3);\n",
       "}\n",
       "\n",
       ".sk-estimator-doc-link.fitted span {\n",
       "  /* fitted */\n",
       "  background: var(--sklearn-color-fitted-level-0);\n",
       "  border: var(--sklearn-color-fitted-level-3);\n",
       "}\n",
       "\n",
       ".sk-estimator-doc-link:hover span {\n",
       "  display: block;\n",
       "}\n",
       "\n",
       "/* \"?\"-specific style due to the `<a>` HTML tag */\n",
       "\n",
       "#sk-container-id-2 a.estimator_doc_link {\n",
       "  float: right;\n",
       "  font-size: 1rem;\n",
       "  line-height: 1em;\n",
       "  font-family: monospace;\n",
       "  background-color: var(--sklearn-color-background);\n",
       "  border-radius: 1rem;\n",
       "  height: 1rem;\n",
       "  width: 1rem;\n",
       "  text-decoration: none;\n",
       "  /* unfitted */\n",
       "  color: var(--sklearn-color-unfitted-level-1);\n",
       "  border: var(--sklearn-color-unfitted-level-1) 1pt solid;\n",
       "}\n",
       "\n",
       "#sk-container-id-2 a.estimator_doc_link.fitted {\n",
       "  /* fitted */\n",
       "  border: var(--sklearn-color-fitted-level-1) 1pt solid;\n",
       "  color: var(--sklearn-color-fitted-level-1);\n",
       "}\n",
       "\n",
       "/* On hover */\n",
       "#sk-container-id-2 a.estimator_doc_link:hover {\n",
       "  /* unfitted */\n",
       "  background-color: var(--sklearn-color-unfitted-level-3);\n",
       "  color: var(--sklearn-color-background);\n",
       "  text-decoration: none;\n",
       "}\n",
       "\n",
       "#sk-container-id-2 a.estimator_doc_link.fitted:hover {\n",
       "  /* fitted */\n",
       "  background-color: var(--sklearn-color-fitted-level-3);\n",
       "}\n",
       "</style><div id=\"sk-container-id-2\" class=\"sk-top-container\"><div class=\"sk-text-repr-fallback\"><pre>add(add(add(add(add(add(Thr1(I39), Thr3(I60)), add(Thr1(I39), Thr2(I60))), add(add(Thr1(I39), Thr3(I60)), add(add(Thr3(I28), Thr1(I39)), add(Thr2(I51), Thr3(I23))))), add(add(add(add(add(add(add(add(Thr3(I39), Thr2(I60)), add(Thr1(I51), Thr3(I23))), add(add(Thr3(I39), Thr1(I60)), add(Thr3(I28), Thr1(I60)))), add(add(add(Thr1(I58), Thr3(I28)), add(Thr2(I28), Thr1(I28))), add(add(0.728, Thr2(I39)), add(Thr1(I51), Thr1(I60))))), add(add(add(add(Thr2(I60), Thr2(I23)), add(Thr2(I51), Thr1(I23))), add(add(Thr1(I28), -0.005), add(Thr1(I28), Thr1(I51)))), add(add(add(Thr1(I60), Thr2(I39)), add(Thr1(I51), Thr3(I60))), add(add(add(add(Thr3(I28), Thr2(I51)), add(Thr2(I51), Thr3(I23))), add(add(Thr3(I39), Thr1(I39)), add(Thr3(I23), Thr1(I39)))), add(add(add(Thr1(I23), Thr3(I28)), add(Thr2(I58), Thr1(I39))), add(add(Thr2(I58), Thr1(I58)), add(Thr1(I51), Thr1(I60)))))))), add(add(Thr2(I51), Thr2(I39)), add(Thr1(I28), Thr3(I28)))), add(Thr2(I58), Thr1(I60))), add(add(Thr1(I39), Thr2(I58)), add(Thr1(I28), Thr3(I23))))), add(add(add(add(Thr1(I39), Thr3(I60)), add(Thr2(I28), Thr1(I51))), add(Thr1(I39), Thr3(I58))), add(add(add(Thr1(I60), Thr1(I28)), add(Thr1(I39), Thr1(I51))), add(add(Thr2(I60), Thr1(I28)), add(Thr1(I23), Thr2(I23)))))), add(add(add(add(add(Thr1(I28), Thr1(I58)), add(Thr2(I23), Thr2(I28))), add(add(Thr1(I39), Thr3(I60)), add(Thr1(I51), Thr3(I39)))), add(add(add(add(add(add(Thr3(I23), Thr1(I28)), add(Thr2(I51), Thr1(I60))), add(add(Thr1(I51), Thr1(I28)), add(Thr3(I60), Thr2(I51)))), add(add(add(0.630, Thr1(I28)), add(add(add(add(add(add(Thr2(I58), Thr1(I58)), add(Thr2(I39), Thr3(I58))), add(add(Thr3(I58), Thr3(I39)), add(Thr2(I58), Thr3(I60)))), add(add(add(Thr1(I51), Thr3(I39)), add(Thr2(I39), Thr1(I60))), add(add(Thr1(I39), Thr3(I28)), add(0.967, 0.434)))), add(add(add(add(Thr1(I28), Thr3(I60)), add(Thr3(I39), Thr1(I39))), add(add(Thr1(I28), Thr3(I51)), add(Thr1(I60), Thr2(I23)))), add(add(add(Thr3(I39), Thr1(I60)), add(Thr3(I28), Thr1(I58))), add(add(Thr3(I28), Thr1(I28)), add(Thr2(I39), Thr1(I39)))))), add(add(add(add(add(-0.797, Thr2(I60)), add(Thr3(I28), Thr1(I51))), add(add(Thr3(I58), Thr2(I28)), add(Thr2(I60), Thr3(I51)))), add(add(add(Thr3(I60), Thr3(I28)), add(Thr1(I28), Thr1(I39))), add(add(Thr1(I51), Thr3(I51)), add(Thr1(I39), Thr3(I51))))), add(add(add(add(Thr1(I58), Thr3(I58)), add(Thr3(I28), Thr1(I39))), add(add(Thr2(I60), Thr1(I28)), add(Thr1(I60), Thr3(I39)))), add(add(add(Thr1(I23), -0.324), add(Thr2(I23), Thr1(I58))), add(add(Thr1(I39), Thr2(I51)), add(Thr1(I60), Thr2(I60)))))))), add(add(Thr3(I39), Thr3(I60)), add(Thr1(I58), Thr1(I28))))), add(Thr3(I28), Thr2(I28))), add(add(Thr2(I51), Thr2(I51)), add(Thr1(I28), Thr1(I58))))), add(add(add(add(Thr1(I28), Thr3(I23)), add(Thr1(I58), Thr3(I58))), add(add(Thr1(I39), Thr3(I58)), add(Thr1(I58), Thr3(I23)))), add(add(add(Thr2(I51), Thr2(I58)), add(Thr2(I23), Thr2(I60))), add(add(Thr1(I58), Thr3(I58)), add(Thr1(I23), Thr2(I39)))))))</pre><b>In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. <br />On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.</b></div><div class=\"sk-container\" hidden><div class=\"sk-item\"><div class=\"sk-estimator fitted sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-2\" type=\"checkbox\" checked><label for=\"sk-estimator-id-2\" class=\"sk-toggleable__label fitted sk-toggleable__label-arrow fitted\">&nbsp;SymbolicRegressor<span class=\"sk-estimator-doc-link fitted\">i<span>Fitted</span></span></label><div class=\"sk-toggleable__content fitted\"><pre>add(add(add(add(add(add(Thr1(I39), Thr3(I60)), add(Thr1(I39), Thr2(I60))), add(add(Thr1(I39), Thr3(I60)), add(add(Thr3(I28), Thr1(I39)), add(Thr2(I51), Thr3(I23))))), add(add(add(add(add(add(add(add(Thr3(I39), Thr2(I60)), add(Thr1(I51), Thr3(I23))), add(add(Thr3(I39), Thr1(I60)), add(Thr3(I28), Thr1(I60)))), add(add(add(Thr1(I58), Thr3(I28)), add(Thr2(I28), Thr1(I28))), add(add(0.728, Thr2(I39)), add(Thr1(I51), Thr1(I60))))), add(add(add(add(Thr2(I60), Thr2(I23)), add(Thr2(I51), Thr1(I23))), add(add(Thr1(I28), -0.005), add(Thr1(I28), Thr1(I51)))), add(add(add(Thr1(I60), Thr2(I39)), add(Thr1(I51), Thr3(I60))), add(add(add(add(Thr3(I28), Thr2(I51)), add(Thr2(I51), Thr3(I23))), add(add(Thr3(I39), Thr1(I39)), add(Thr3(I23), Thr1(I39)))), add(add(add(Thr1(I23), Thr3(I28)), add(Thr2(I58), Thr1(I39))), add(add(Thr2(I58), Thr1(I58)), add(Thr1(I51), Thr1(I60)))))))), add(add(Thr2(I51), Thr2(I39)), add(Thr1(I28), Thr3(I28)))), add(Thr2(I58), Thr1(I60))), add(add(Thr1(I39), Thr2(I58)), add(Thr1(I28), Thr3(I23))))), add(add(add(add(Thr1(I39), Thr3(I60)), add(Thr2(I28), Thr1(I51))), add(Thr1(I39), Thr3(I58))), add(add(add(Thr1(I60), Thr1(I28)), add(Thr1(I39), Thr1(I51))), add(add(Thr2(I60), Thr1(I28)), add(Thr1(I23), Thr2(I23)))))), add(add(add(add(add(Thr1(I28), Thr1(I58)), add(Thr2(I23), Thr2(I28))), add(add(Thr1(I39), Thr3(I60)), add(Thr1(I51), Thr3(I39)))), add(add(add(add(add(add(Thr3(I23), Thr1(I28)), add(Thr2(I51), Thr1(I60))), add(add(Thr1(I51), Thr1(I28)), add(Thr3(I60), Thr2(I51)))), add(add(add(0.630, Thr1(I28)), add(add(add(add(add(add(Thr2(I58), Thr1(I58)), add(Thr2(I39), Thr3(I58))), add(add(Thr3(I58), Thr3(I39)), add(Thr2(I58), Thr3(I60)))), add(add(add(Thr1(I51), Thr3(I39)), add(Thr2(I39), Thr1(I60))), add(add(Thr1(I39), Thr3(I28)), add(0.967, 0.434)))), add(add(add(add(Thr1(I28), Thr3(I60)), add(Thr3(I39), Thr1(I39))), add(add(Thr1(I28), Thr3(I51)), add(Thr1(I60), Thr2(I23)))), add(add(add(Thr3(I39), Thr1(I60)), add(Thr3(I28), Thr1(I58))), add(add(Thr3(I28), Thr1(I28)), add(Thr2(I39), Thr1(I39)))))), add(add(add(add(add(-0.797, Thr2(I60)), add(Thr3(I28), Thr1(I51))), add(add(Thr3(I58), Thr2(I28)), add(Thr2(I60), Thr3(I51)))), add(add(add(Thr3(I60), Thr3(I28)), add(Thr1(I28), Thr1(I39))), add(add(Thr1(I51), Thr3(I51)), add(Thr1(I39), Thr3(I51))))), add(add(add(add(Thr1(I58), Thr3(I58)), add(Thr3(I28), Thr1(I39))), add(add(Thr2(I60), Thr1(I28)), add(Thr1(I60), Thr3(I39)))), add(add(add(Thr1(I23), -0.324), add(Thr2(I23), Thr1(I58))), add(add(Thr1(I39), Thr2(I51)), add(Thr1(I60), Thr2(I60)))))))), add(add(Thr3(I39), Thr3(I60)), add(Thr1(I58), Thr1(I28))))), add(Thr3(I28), Thr2(I28))), add(add(Thr2(I51), Thr2(I51)), add(Thr1(I28), Thr1(I58))))), add(add(add(add(Thr1(I28), Thr3(I23)), add(Thr1(I58), Thr3(I58))), add(add(Thr1(I39), Thr3(I58)), add(Thr1(I58), Thr3(I23)))), add(add(add(Thr2(I51), Thr2(I58)), add(Thr2(I23), Thr2(I60))), add(add(Thr1(I58), Thr3(I58)), add(Thr1(I23), Thr2(I39)))))))</pre></div> </div></div></div></div>"
      ],
      "text/plain": [
       "SymbolicRegressor(feature_names=Index(['Thr1(I23)', 'Thr2(I23)', 'Thr3(I23)', 'Thr1(I28)', 'Thr2(I28)',\n",
       "       'Thr3(I28)', 'Thr1(I39)', 'Thr2(I39)', 'Thr3(I39)', 'Thr1(I51)',\n",
       "       'Thr2(I51)', 'Thr3(I51)', 'Thr1(I58)', 'Thr2(I58)', 'Thr3(I58)',\n",
       "       'Thr1(I60)', 'Thr2(I60)', 'Thr3(I60)'],\n",
       "      dtype='object'),\n",
       "                  function_set=['add'], max_samples=0.9, p_crossover=0.7,\n",
       "                  p_hoist_mutation=0.03132735169322752,\n",
       "                  p_point_mutation=0.052454815957287154,\n",
       "                  p_subtree_mutation=0.08463109166860174,\n",
       "                  parsimony_coefficient=0.01, population_size=10000,\n",
       "                  random_state=10, stopping_criteria=0.01, verbose=1)"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function_set = ['add']\n",
    "\n",
    "# Run optimization based on the task\n",
    "if task == 0:\n",
    "    result = gp_minimize(sym_regression, pbounds_regression, n_calls=20, random_state=1, verbose=True)\n",
    "    \n",
    "    # Print and use the best parameters for Task 0\n",
    "    print(\"Best parameters for Task 0:\", result.x)\n",
    "    best_params = result.x\n",
    "    est_gp = SR(\n",
    "        population_size=10000,\n",
    "        function_set=function_set,\n",
    "        generations=20,\n",
    "        stopping_criteria=0.01,\n",
    "        p_crossover=0.7,\n",
    "        p_subtree_mutation=best_params[0],\n",
    "        p_hoist_mutation=best_params[1],\n",
    "        p_point_mutation=best_params[2],\n",
    "        max_samples=0.9,\n",
    "        verbose=1,\n",
    "        metric=\"mean absolute error\",\n",
    "        parsimony_coefficient=0.01,\n",
    "        random_state=10,\n",
    "        feature_names=x_train.columns\n",
    "    )\n",
    "\n",
    "elif task == 1:\n",
    "    result = gp_minimize(sym_classification, pbounds_classification, n_calls=20, random_state=1, verbose=True)\n",
    "    \n",
    "    # Print and use the best parameters for Task 1\n",
    "    print(\"Best parameters for Task 1:\", result.x)\n",
    "    best_params = result.x\n",
    "    est_gp = SC(\n",
    "        population_size=20000,\n",
    "        function_set=function_set,\n",
    "        generations=100,\n",
    "        stopping_criteria=10**(-4),\n",
    "        p_crossover=0.7,\n",
    "        p_subtree_mutation=best_params[0],\n",
    "        p_hoist_mutation=best_params[1],\n",
    "        p_point_mutation=best_params[2],\n",
    "        max_samples=0.9,\n",
    "        verbose=1,\n",
    "        random_state=1,\n",
    "        feature_names=x_train.columns\n",
    "    )            \n",
    "\n",
    "# Fit the model with the optimized parameters\n",
    "est_gp.fit(x_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79a23c05",
   "metadata": {},
   "source": [
    "# Check for overfitting: comparing performance between training and test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "c777f872",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training MAE: 10.288735017546117\n",
      "Test MAE: 10.635138766312364\n"
     ]
    }
   ],
   "source": [
    "# Predict on training and test sets\n",
    "if task == 0:\n",
    "    train_predictions = est_gp.predict(x_train)\n",
    "    test_predictions = est_gp.predict(x_test)\n",
    "\n",
    "    train_mae = mean_absolute_error(y_train, train_predictions)\n",
    "    test_mae = mean_absolute_error(y_test, test_predictions)\n",
    "\n",
    "    print(f\"Training MAE: {train_mae}\")\n",
    "    print(f\"Test MAE: {test_mae}\")\n",
    "    \n",
    "elif task == 1:\n",
    "    train_predictions = est_gp.predict_proba(x_train)\n",
    "    test_predictions = est_gp.predict_proba(x_test)\n",
    "    fpr_sym_train, tpr_sym_train, thresholds_train = roc_curve(y_train, train_predictions[:, 1])\n",
    "    roc_auc_sym_train = auc(fpr_sym_train, tpr_sym_train)\n",
    "    fpr_sym_test, tpr_sym_test, thresholds_test = roc_curve(y_test, test_predictions[:, 1])\n",
    "    roc_auc_sym_test = auc(fpr_sym_test, tpr_sym_test)\n",
    "    \n",
    "    print(f\"Training AUC: {roc_auc_sym_train}\")\n",
    "    print(f\"Test AUC: {roc_auc_sym_test}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "e54f1237",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The model is generalizing well.\n"
     ]
    }
   ],
   "source": [
    "# Threshold ratio for overfitting\n",
    "if task==0:\n",
    "    overfit_ratio = 1.5 \n",
    "\n",
    "    # Check for overfitting\n",
    "    if test_mae > overfit_ratio * train_mae:\n",
    "        print(\"Warning: The model might be overfitting!\")\n",
    "    else:\n",
    "        print(\"The model is generalizing well.\")\n",
    "\n",
    "elif task ==1:\n",
    "    # Set a threshold for AUC difference\n",
    "    auc_diff_threshold = 0.1\n",
    "\n",
    "    # Calculate the absolute difference between train and test AUC\n",
    "    auc_diff = abs(roc_auc_sym_train - roc_auc_sym_test)\n",
    "\n",
    "    # Issue a warning if the difference exceeds the threshold\n",
    "    if auc_diff > auc_diff_threshold:\n",
    "        print(f\"Warning: Significant difference in AUC between training and test sets (Difference: {auc_diff:.3f}). The model might be overfitting!\")\n",
    "    else:\n",
    "        print(\"The model is generalizing well.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6bbd72c2",
   "metadata": {},
   "source": [
    "# Check for performance on different subsets: Cross-validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32e117d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "if task==0:\n",
    "    cv = KFold(n_splits=5, shuffle=True, random_state=42)\n",
    "\n",
    "    # Perform cross-validation with the defined cross-validator\n",
    "    scores = cross_val_score(est_gp, X, y, cv=cv, scoring='neg_mean_absolute_error')\n",
    "    mean_cv_score = -np.mean(scores)\n",
    "    std_cv_score = np.std(scores)\n",
    "    \n",
    "    # Print cross-validated AUC result\n",
    "    print(f\"Cross-validated MAE: {mean_cv_score}\")\n",
    "    print(f\"Standard deviation of Cross-validated MAE: {std_cv_score}\")\n",
    "\n",
    "elif task ==1:\n",
    "    cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)\n",
    "\n",
    "    # Get cross-validated predictions (probabilities) on the entire dataset\n",
    "    cv_predictions = cross_val_predict(est_gp, X, y, cv=cv, method='predict_proba')\n",
    "\n",
    "    # Calculate ROC AUC for the cross-validated predictions\n",
    "    fpr_cv, tpr_cv, _ = roc_curve(y, cv_predictions[:, 1])\n",
    "    roc_auc_cv = auc(fpr_cv, tpr_cv)\n",
    "\n",
    "    # Print cross-validated AUC result\n",
    "    print(f\"Cross-validated AUC: {roc_auc_cv}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00f05aa0",
   "metadata": {},
   "outputs": [],
   "source": [
    "if task ==0:\n",
    "    print(f\"Cross-validated scores: {-scores}\")\n",
    "    print(f\"Mean cross-validated score: {mean_cv_score}\")\n",
    "    print(f\"Standard deviation of CV scores: {std_cv_score}\")\n",
    "\n",
    "    train_predictions = est_gp.predict(x_train)\n",
    "    train_mae = mean_absolute_error(y_train, train_predictions)\n",
    "\n",
    "    print(f\"Training MAE: {train_mae}\")\n",
    "\n",
    "    overfit_threshold = 1.5  # You can adjust this value\n",
    "\n",
    "    # Check for overfitting\n",
    "    if mean_cv_score > overfit_threshold * train_mae:\n",
    "        print(\"Warning: The model might be overfitting! The cross-validation error is significantly higher than the training error.\")\n",
    "    elif std_cv_score > 0.15 * mean_cv_score:  # Check for high variance in CV scores \n",
    "        print(\"Warning: The model might be unstable! High variance in cross-validation performance.\")\n",
    "    else:\n",
    "        print(\"The model is generalizing well.\")\n",
    "    \n",
    "elif task ==1:\n",
    "    # Calculate AUROC for each fold\n",
    "    cv_aucs = cross_val_score(est_gp, X, y, cv=cv, scoring='roc_auc')\n",
    "\n",
    "    # Print the AUROC for each fold\n",
    "    print(f\"AUROC for each fold: {cv_aucs}\")\n",
    "\n",
    "    # Calculate the mean and standard deviation of AUROC across folds\n",
    "    mean_auc = np.mean(cv_aucs)\n",
    "    std_auc = np.std(cv_aucs)\n",
    "\n",
    "    # Print average and standard deviation\n",
    "    print(f\"Mean AUROC: {mean_auc}\")\n",
    "    print(f\"Standard deviation of AUROC: {std_auc}\")\n",
    "\n",
    "    # Set thresholds for acceptable performance\n",
    "    min_auc_threshold = 0.7  # Minimum acceptable AUROC\n",
    "    max_std_threshold = 0.1  # Maximum acceptable standard deviation of AUROC across folds\n",
    "\n",
    "    # Check for poor performance\n",
    "    if mean_auc < min_auc_threshold:\n",
    "        print(f\"Warning: Low mean AUROC ({mean_auc:.3f}) across folds, indicating poor model performance.\")\n",
    "\n",
    "    if std_auc > max_std_threshold:\n",
    "        print(f\"Warning: High variability in AUROC across folds (Standard Deviation: {std_auc:.3f}), indicating model instability.\")\n",
    "    else:\n",
    "        print(\"The model is generalizing well.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34927a90",
   "metadata": {},
   "source": [
    "# Define a Python dictionary \"converter\"\n",
    "\n",
    "The converter dictionary provides a mapping between string keys (e.g., 'sub', 'div', 'mul', 'add') and lambda functions representing arithmetic operations. This dictionary can be used to select the appropriate operation based on a given key and apply it to numerical values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "c0a8f45e",
   "metadata": {},
   "outputs": [],
   "source": [
    "converter = {\n",
    "    'sub': lambda x, y : x - y,\n",
    "    'div': lambda x, y : x/y,\n",
    "    'mul': lambda x, y : x*y,\n",
    "    'add': lambda x, y : x + y,\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5110a7f7",
   "metadata": {},
   "source": [
    "# Evaluate the performance of your symbolic regression model (est_gp)\n",
    "\n",
    "est_gp.score(x_test, y_test): calculates the R^2 score of the symbolic regression model (est_gp) on the test data (x_test, y_test). The R^2 score measures the proportion of the variance in the target variable that is predictable from the features. Higher R^2 scores indicate a better fit of the model to the data.\n",
    "\n",
    "next_e: converts the best individual (mathematical expression) learned by the symbolic regression model (est_gp._program) into a SymPy symbolic expression. \n",
    "\n",
    "    sympify is a function from SymPy that converts a string representing a mathematical expression into a SymPy\n",
    "    expression. The locals parameter is used to provide the converter dictionary, which maps mathematical\n",
    "    operations to corresponding lambda functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "91fe2a1b",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R2 score: 0.7375566846336077\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 16 \\operatorname{Thr}_{1}{\\left(I_{23} \\right)} + 13 \\operatorname{Thr}_{1}{\\left(I_{28} \\right)} + 15 \\operatorname{Thr}_{1}{\\left(I_{39} \\right)} + 8 \\operatorname{Thr}_{1}{\\left(I_{51} \\right)} + 13 \\operatorname{Thr}_{1}{\\left(I_{58} \\right)} + 17 \\operatorname{Thr}_{1}{\\left(I_{60} \\right)} + 6 \\operatorname{Thr}_{2}{\\left(I_{23} \\right)} + 5 \\operatorname{Thr}_{2}{\\left(I_{28} \\right)} + 7 \\operatorname{Thr}_{2}{\\left(I_{39} \\right)} + 5 \\operatorname{Thr}_{2}{\\left(I_{51} \\right)} + 13 \\operatorname{Thr}_{2}{\\left(I_{58} \\right)} + 8 \\operatorname{Thr}_{2}{\\left(I_{60} \\right)} + 8 \\operatorname{Thr}_{3}{\\left(I_{23} \\right)} + 10 \\operatorname{Thr}_{3}{\\left(I_{28} \\right)} + 5 \\operatorname{Thr}_{3}{\\left(I_{39} \\right)} + 7 \\operatorname{Thr}_{3}{\\left(I_{51} \\right)} + 8 \\operatorname{Thr}_{3}{\\left(I_{58} \\right)} + 5 \\operatorname{Thr}_{3}{\\left(I_{60} \\right)} + 10 \\operatorname{Thr}_{4}{\\left(I_{23} \\right)} + 12 \\operatorname{Thr}_{4}{\\left(I_{28} \\right)} + 7 \\operatorname{Thr}_{4}{\\left(I_{39} \\right)} + 8 \\operatorname{Thr}_{4}{\\left(I_{51} \\right)} + 13 \\operatorname{Thr}_{4}{\\left(I_{58} \\right)} + 10 \\operatorname{Thr}_{4}{\\left(I_{60} \\right)} + 3.206$"
      ],
      "text/plain": [
       "16*Thr1(I23) + 13*Thr1(I28) + 15*Thr1(I39) + 8*Thr1(I51) + 13*Thr1(I58) + 17*Thr1(I60) + 6*Thr2(I23) + 5*Thr2(I28) + 7*Thr2(I39) + 5*Thr2(I51) + 13*Thr2(I58) + 8*Thr2(I60) + 8*Thr3(I23) + 10*Thr3(I28) + 5*Thr3(I39) + 7*Thr3(I51) + 8*Thr3(I58) + 5*Thr3(I60) + 10*Thr4(I23) + 12*Thr4(I28) + 7*Thr4(I39) + 8*Thr4(I51) + 13*Thr4(I58) + 10*Thr4(I60) + 3.206"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print('R2 score:',est_gp.score(x_test,y_test))\n",
    "next_e = sympify(str(est_gp._program), locals = converter)\n",
    "next_e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "896b2d71",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Process the symbolic expression (next_e) to extract and analyze its terms, coefficients, and values \n",
    "# associated with specific variables. \n",
    "# After processing, dataframes are created to store the results and then save them to Excel files.\n",
    "\n",
    "if task ==0:\n",
    "    e_str = str(next_e).replace(\" \", \"\")\n",
    "    terms = e_str.split(\"+\")\n",
    "    variables = selected_categorical_variables\n",
    "    term_grouped = {str(var): [] for var in variables}\n",
    "    term_coefs = {str(var): [0]*len(categorical_response_list) for var in variables}\n",
    "    term_val = {f\"Q{var}\": [0]*len(categorical_response_list) for var in variables}\n",
    "        \n",
    "        \n",
    "    for term in terms:\n",
    "        if \"I\" not in term:\n",
    "            continue\n",
    "        left = term.index('I')\n",
    "        for t in range(len(variables)):\n",
    "            if term[left+1:-1] == str(variables[t]):\n",
    "                term_grouped[str(variables[t])].append(term)\n",
    "\n",
    "    \n",
    "\n",
    "    for key in term_grouped.keys():\n",
    "        lst = term_grouped[key]\n",
    "        for item in lst:\n",
    "            idx = int(item[item.index(\"r\")+1:item.index(\"(\")]) - 1\n",
    "            coef = 1\n",
    "            if item[0]!=\"T\":\n",
    "                coef = int(item[:item.index(\"*\")])\n",
    "            term_coefs[key][idx] = coef\n",
    "\n",
    "        for d in range(len(categorical_response_list)):\n",
    "            term_val['Q'+key][d] = term_val['Q'+key][d-1] + term_coefs[key][d]\n",
    "\n",
    "    \n",
    "elif task==1:\n",
    "    next_e_str = str(next_e).replace(\" \", \"\")\n",
    "    terms = re.split('\\+|-', next_e_str)\n",
    "\n",
    "    term_grouped = {}\n",
    "    term_coefs = {}\n",
    "    \n",
    "\n",
    "    for term in terms:\n",
    "        if '(' in term and ')' in term:\n",
    "            start_idx = term.index('(') + 1\n",
    "            end_idx = term.index(')')\n",
    "            var = term[start_idx:end_idx]\n",
    "        else:\n",
    "            if '*' in term:\n",
    "                var = term.split('*')[-1]\n",
    "            else:\n",
    "                continue\n",
    "        \n",
    "        if var not in term_grouped:\n",
    "            term_grouped[var] = []\n",
    "        term_grouped[var].append(term)\n",
    "       \n",
    "        if '*' in term:\n",
    "            coef = int(term.split('*')[0])\n",
    "        else:\n",
    "            coef = 1\n",
    "\n",
    "        if var not in term_coefs:\n",
    "            term_coefs[var] = []\n",
    "        term_coefs[var].append(coef)\n",
    "\n",
    "    max_numb = {var: 0 for var in term_grouped.keys()}\n",
    "\n",
    "    for var, terms_list in term_grouped.items():\n",
    "        for term in terms_list:\n",
    "            num_match = re.search(r'\\d+(?=\\()', term)\n",
    "            if num_match:\n",
    "                num = int(num_match.group(0))\n",
    "                if num > max_numb[var]:\n",
    "                    max_numb[var] = num\n",
    "\n",
    "    maximum_row = int(max(max_numb.values())) + 1\n",
    "    term_val = {key: [0]*maximum_row for key in term_coefs.keys()}\n",
    "\n",
    "    for var, terms_list in term_grouped.items():\n",
    "        for term in terms_list:\n",
    "            if '(' not in term:\n",
    "                if '*' in term:\n",
    "                    coef = int(term.split('*')[0])\n",
    "                    for i in range(maximum_row):\n",
    "                        term_val[var][i] = coef if var in term else 0\n",
    "                else:\n",
    "                    for i in range(maximum_row):\n",
    "                        term_val[var][i] = 1 if var in term else 0\n",
    "            else:\n",
    "                num_match = re.search(r'\\d+(?=\\()', term)\n",
    "                if num_match:\n",
    "                    num = int(num_match.group(0))\n",
    "                else:\n",
    "                    num = 1\n",
    "            \n",
    "                coef = 1\n",
    "                if '*' in term:\n",
    "                    coef = int(term.split('*')[0])\n",
    "                \n",
    "                if 'iThr' in term:\n",
    "                    for i in range(0, min(num, maximum_row)):\n",
    "                        term_val[var][i] += coef\n",
    "\n",
    "                elif 'Thr' in term:\n",
    "                    for i in range(num, maximum_row):\n",
    "                        term_val[var][i] += coef\n",
    "    \n",
    "    table0 = pd.DataFrame(term_val)\n",
    "    print(table0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "9e13b703",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   Q23  Q28  Q39  Q51  Q58  Q60\n",
      "0   16   13   15    8   13   17\n",
      "1   22   18   22   13   26   25\n",
      "2   30   28   27   20   34   30\n",
      "3   40   40   34   28   47   40\n",
      "   23  28  39  51  58  60\n",
      "0  16  13  15   8  13  17\n",
      "1   6   5   7   5  13   8\n",
      "2   8  10   5   7   8   5\n",
      "3  10  12   7   8  13  10\n",
      "229\n"
     ]
    }
   ],
   "source": [
    "if task == 0:\n",
    "    table0 = pd.DataFrame(term_val)\n",
    "    table0.to_excel('SCORES/table_temp.xlsx')\n",
    "\n",
    "    wei = pd.DataFrame(term_coefs)\n",
    "    wei.to_excel('SCORES/partial_weights.xlsx')\n",
    "\n",
    "    print(table0)\n",
    "    print(wei)\n",
    "    print(table0.iloc[-1].sum())\n",
    "    \n",
    "elif task ==1:\n",
    "    results = {}\n",
    "    max_values = {}\n",
    "\n",
    "    \n",
    "    for column in table0.columns:\n",
    "        \n",
    "        column_results = {}\n",
    "        unique_values = table0[column].unique()\n",
    "\n",
    "        if pd.api.types.is_numeric_dtype(table0[column]):\n",
    "            max_values[column] = table0[column].max()\n",
    "\n",
    "        for unique_value in unique_values:\n",
    "            indices = []\n",
    "\n",
    "            for i in range(len(table0)):\n",
    "                if table0.iloc[i][column] == unique_value:\n",
    "                    indices.append(i)\n",
    "\n",
    "            if indices: \n",
    "                max_index = max(indices)\n",
    "                min_index = min(indices)\n",
    "                column_results[unique_value] = (min_index, max_index)\n",
    "            else:\n",
    "                column_results[unique_value] = (None, None)\n",
    "\n",
    "        results[column] = column_results\n",
    "\n",
    "    all_scores = []\n",
    "\n",
    "    for column, column_result in results.items():\n",
    "        print(f\"Variable: {column}\")\n",
    "        actual_max_value = max_values.get(column, None)\n",
    "        scores = []\n",
    "\n",
    "        if column in binary_variables:\n",
    "            unique_value = list(column_result.keys())[0]\n",
    "            print(f\"Score: {unique_value}\")\n",
    "            scores.append([column, \"\", unique_value])\n",
    "        else:\n",
    "            last_range_start = None\n",
    "            last_range_value = None\n",
    "\n",
    "            for value, (min_index, max_index) in column_result.items():\n",
    "                if min_index is not None and max_index is not None:\n",
    "                    if value == actual_max_value:\n",
    "                        last_range_start = min_index\n",
    "                        last_range_value = value\n",
    "                        break\n",
    "\n",
    "            for value, (min_index, max_index) in column_result.items():\n",
    "                if min_index is not None and max_index is not None:\n",
    "                    if value == last_range_value:\n",
    "                        range_str = f\"{min_index}-{max_index}\"\n",
    "                    elif max_index == maximum_row-1:\n",
    "                        range_str = f\"{min_index}-{max(df[column])}\"  \n",
    "                    else:\n",
    "                        range_str = f\"{min_index}-{max_index}\"\n",
    "                        \n",
    "                    print(f\"Range: {range_str}   Score: {value}\")\n",
    "\n",
    "                    scores.append([column, range_str, value])\n",
    "\n",
    "        print()  \n",
    "\n",
    "        all_scores.extend(scores)\n",
    "\n",
    "    final_score_table = pd.DataFrame(all_scores, columns=['Variable', 'Range', 'Score'])\n",
    "\n",
    "    final_score_table['Range'] = final_score_table.apply(lambda row: re.sub(r'\\d+-$','4-', row['Range']) if row['Variable'].startswith('ISI') else row['Range'], axis=1)\n",
    "\n",
    "    final_score_table.to_excel('SCORES/score_table_classification.xlsx', sheet_name='Results', index=False)\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "0b35ebaf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mean absolute error\n",
    "def mae(y_true, predictions):\n",
    "    y_true, predictions = np.array(y_true), np.array(predictions)\n",
    "    return np.mean(np.abs(y_true - predictions))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "0551327a",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "db898da183a945d29e367f37af58a634",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "진행률:   0%|          | 0/24 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10.010752688172044\n",
      "   Q23  Q28  Q39  Q51  Q58  Q60\n",
      "0   16   13   15    8   13   17\n",
      "1   22   18   22   13   26   25\n",
      "2   30   28   27   20   34   30\n",
      "3   40   40   34   28   47   51\n"
     ]
    }
   ],
   "source": [
    "if task ==0:\n",
    "    # Prepare data\n",
    "    X2 = df_linear[:][1:]\n",
    "    y2 = df[categorical_items][:][1:].sum(axis=1).values\n",
    "    x_train2, x_test2, y_train2, y_test2 = train_test_split(X2, y2, test_size=0.3, random_state=10)\n",
    "\n",
    "    # Deep copy the test data\n",
    "    xtest = copy.deepcopy(x_test2)\n",
    "\n",
    "    # Load weights and cumulative term values\n",
    "    wei = pd.read_excel('SCORES/partial_weights.xlsx')\n",
    "    table0 = pd.read_excel('SCORES/table_temp.xlsx')\n",
    "\n",
    "    # Calculate difference between target value and sum of last row in table0\n",
    "    candy = max_items*len(categorical_response_list) - table0[table0.columns[1:]].iloc[-1].sum()\n",
    "\n",
    "    # Fine-tune model if difference is not zero\n",
    "    if candy != 0:\n",
    "        # Generate cases to explore different adjustments to the coefficients\n",
    "        cases = list(itertools.product(table0.columns[1:], [0, 1, 2, 3]))\n",
    "        table_opt = 0\n",
    "        MAE_opt = 9999999\n",
    "\n",
    "        # Iterate through cases\n",
    "        for case in tqdm(cases, desc='진행률'):\n",
    "            # Deep copy weights\n",
    "            wgt = copy.deepcopy(wei)\n",
    "\n",
    "            col, row = case[0], case[1]\n",
    "            wgt[col[1:]][row] += candy\n",
    "\n",
    "            # Skip if weight becomes negative\n",
    "            if wgt[col[1:]][row] < 0:\n",
    "                continue\n",
    "\n",
    "            # Calculate cumulative term values with adjusted weights\n",
    "    \n",
    "            term_val = {f\"Q{var}\": [0]*len(categorical_response_list) for var in selected_categorical_variables}\n",
    "            for key in wei.columns[1:]:\n",
    "                for d in range(4):\n",
    "                    term_val['Q' + key][d] = term_val['Q' + key][d - 1] + wgt[key][d]\n",
    "\n",
    "            weight = pd.DataFrame(term_val)\n",
    "\n",
    "            # Apply adjusted weights to test data\n",
    "            dictlist = []\n",
    "            for key in selected_categorical_variables:\n",
    "                dic = (key, {k: weight[\"Q\" + str(key)].iloc[k] for k in range(len(categorical_response_list))})\n",
    "                dictlist.append(dic)\n",
    "            ttsdic = dict(dictlist)\n",
    "\n",
    "            for key in selected_categorical_variables:\n",
    "                xtest[key] = x_test2[key].apply(lambda x: ttsdic[key][x - categorical_response_list[0]])\n",
    "\n",
    "            # Calculate mean absolute error\n",
    "            y_pred2 = xtest.sum(axis=1)\n",
    "            MAE = mae(y_test2, y_pred2)\n",
    "\n",
    "            # Update optimal solution if MAE is lower\n",
    "            if MAE < MAE_opt:\n",
    "                MAE_opt = MAE\n",
    "                table_opt = copy.deepcopy(weight)\n",
    "            time.sleep(0.0001)\n",
    "\n",
    "        # Print and save results\n",
    "        print(MAE_opt)\n",
    "        table_opt.to_excel('SCORES/Scores_regression_FineTuned.xlsx')\n",
    "        print(table_opt)\n",
    "        \n",
    "elif task ==1:\n",
    "    thres_unprocessed = float(terms[-1])\n",
    "    thres = thres_unprocessed*len(all_items)\n",
    "    print(f\"If total score is greater than {thres:.2f}, then the individual has the disease/disorder.\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ea959ff-446e-4029-9f11-dc00845711b2",
   "metadata": {},
   "source": [
    "# Saving the Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3472e8d9-33b3-42b7-8519-5a0933183e5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('sym_reg_model_mcqi.pkl', 'wb') as file:\n",
    "    pickle.dump(est_gp, file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17168199-5cb1-48bf-a84e-7bc7ff858779",
   "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}