{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Genetic algorithm to optimize thresholds + analysis\n\nSince space of thresholds is very vast, techniques like random search won't work to optimize overall f1.\n\nThis notebook presents genetic algorithm to apply small threshold shifts to optimize overall f1.","metadata":{}},{"cell_type":"code","source":"import os\nimport pickle\nimport math\n\nfrom multiprocessing import Pool\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom sklearn.metrics import roc_curve, auc, f1_score, mean_squared_error","metadata":{"execution":{"iopub.status.busy":"2023-06-02T11:41:22.192943Z","iopub.execute_input":"2023-06-02T11:41:22.193292Z","iopub.status.idle":"2023-06-02T11:41:22.198690Z","shell.execute_reply.started":"2023-06-02T11:41:22.193263Z","shell.execute_reply":"2023-06-02T11:41:22.197463Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with open('/kaggle/input/obtained-data/statistics.pkl', 'rb') as f:\n    true, model_predictions = pickle.load(f)","metadata":{"execution":{"iopub.status.busy":"2023-06-02T11:41:22.201318Z","iopub.execute_input":"2023-06-02T11:41:22.201849Z","iopub.status.idle":"2023-06-02T11:41:22.219926Z","shell.execute_reply.started":"2023-06-02T11:41:22.201812Z","shell.execute_reply":"2023-06-02T11:41:22.219125Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Shared threshol\n- lets find single threshold for all questions, that yields best score","metadata":{}},{"cell_type":"code","source":"scores = []\nthresholds = []\nbest_score = 0\nbest_threshold_global = 0\n\nfor threshold in np.arange(0, 1, 0.01):\n    preds = (model_predictions.values.reshape(-1) > threshold).astype('int')\n    m = f1_score(true.values.reshape(-1), preds, average='macro')\n    scores.append(m)\n    thresholds.append(threshold)\n    if m > best_score:\n        best_score = m\n        best_threshold_global = threshold\n\nprint('best score: ', best_score)\nprint('best threshold: ', best_threshold_global)","metadata":{"execution":{"iopub.status.busy":"2023-06-02T11:41:22.221156Z","iopub.execute_input":"2023-06-02T11:41:22.221560Z","iopub.status.idle":"2023-06-02T11:41:32.946125Z","shell.execute_reply.started":"2023-06-02T11:41:22.221537Z","shell.execute_reply":"2023-06-02T11:41:32.945019Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# plot the f1 score\nplt.plot(thresholds, scores)\nplt.xlabel('threshold')\nplt.ylabel('f1 score')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-02T11:41:32.947348Z","iopub.execute_input":"2023-06-02T11:41:32.947789Z","iopub.status.idle":"2023-06-02T11:41:33.094616Z","shell.execute_reply.started":"2023-06-02T11:41:32.947746Z","shell.execute_reply":"2023-06-02T11:41:33.093900Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# data prep\ndisplay(model_predictions.head())\ndisplay(true.head())\n\npred = []\nlabels = []\nfor i in range(18):\n    pred.append(np.array(model_predictions[i].values.reshape(-1)))\n    print(pred[i])\n    labels.append(np.array(true[i].values.reshape(-1)))\n    print(labels[i])\n","metadata":{"execution":{"iopub.status.busy":"2023-06-02T11:41:33.097746Z","iopub.execute_input":"2023-06-02T11:41:33.099286Z","iopub.status.idle":"2023-06-02T11:41:33.138754Z","shell.execute_reply.started":"2023-06-02T11:41:33.099258Z","shell.execute_reply":"2023-06-02T11:41:33.137719Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Separated thresholds\n- let's analyze thresholds for each question individually","metadata":{}},{"cell_type":"code","source":"thresholds_separated = []\nfor i, (pred, lab) in enumerate(zip(pred, labels)):\n    print(f'question {i+1}')\n    # plot the roc curve\n    fpr, tpr, threshold = roc_curve(lab, pred)\n    # thresholds_separated.append(threshold)\n    roc_auc = auc(fpr, tpr)\n    plt.plot(fpr, tpr, label=f'question {i+1} (AUC = {roc_auc:.2f})')\n    plt.plot([0, 1], [0, 1], 'b--')\n    plt.legend()\n    plt.xlabel('False Positive Rate')\n    plt.ylabel('True Positive Rate')\n    plt.show()\n    print(f'AUC: {roc_auc:.2f}')\n    # calculate optimal threshold based on f1 score\n    scores = []\n    thresholds = []\n    best_score = 0\n    best_threshold = 0\n    for threshold in np.arange(0, 1, 0.01):\n        preds = (pred > threshold).astype('int')\n        m = f1_score(lab, preds, average='macro')\n        scores.append(m)\n        thresholds.append(threshold)\n        if m > best_score:\n            best_score = m\n            best_threshold = threshold\n    # plot the f1 score\n    plt.plot(thresholds, scores)\n    plt.xlabel('threshold')\n    plt.ylabel('f1 score')\n    plt.show()\n    print('best score: ', best_score)\n    print('best threshold: ', best_threshold)\n    thresholds_separated.append(best_threshold)\nprint(thresholds_separated)","metadata":{"execution":{"iopub.status.busy":"2023-06-02T11:41:33.140587Z","iopub.execute_input":"2023-06-02T11:41:33.140892Z","iopub.status.idle":"2023-06-02T11:41:49.073021Z","shell.execute_reply.started":"2023-06-02T11:41:33.140867Z","shell.execute_reply":"2023-06-02T11:41:49.072157Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Genetic algorithm\nfinally implementation of genetic algorithm\n- tournament selection since fitness of all individuals is very similar for roulette selection\n- rate = how often individual is changed\n- impact = how much individual is changed","metadata":{}},{"cell_type":"code","source":"best_scores = []\noutput_file = 'thresholds.txt'\n# Population size and other parameters\npopulation_size = 500\nnum_generations = 1000\nearly_stopping = 20\nelite_size = population_size // 20\ntournament_size = population_size // 50\nmutation_rate = 0.8\nmutation_impact = 0.1\ncrossover_rate = 0.2\ncrossover_impact = 0.1\n\nprint(f'Population size: {population_size}')\nprint(f'Number of generations: {num_generations}')\nprint(f'Early stopping: {early_stopping}')\nprint(f'Elite size: {elite_size}')\nprint(f'Tournament size: {tournament_size}')\nprint()\n\nwith open(output_file, 'w') as f:\n    f.write('score,rmse,thresholds\\n')\n\n# Initialize population\noriginal_thresholds = [best_threshold_global for _ in range(18)]\npopulation = [original_thresholds.copy() for _ in range(population_size)]\n\n# Evaluate fitness of individuals\ndef evaluate_fitness(individual):\n    \"\"\"evaluate fitness of an individual - f1 score on predictions\"\"\"\n    preds = model_predictions.copy()\n    for i in range(18):\n        preds[i] = (model_predictions[i] > individual[i]).astype('int')\n    return f1_score(true.values.reshape(-1), preds.values.reshape(-1), average='macro')\n\n\n# Averaging crossover\ndef average_crossover(parent1, parent2):\n    \"\"\"crossover that averages random genes from the two parents\"\"\"\n    child = parent1.copy()\n    for i, (p1, p2) in enumerate(zip(parent1, parent2)):\n        if np.random.random() < crossover_impact:\n            child[i] = (p1 + p2) / 2\n    return child\n\n\n# Selective crossover\ndef selective_crossover(parent1, parent2):\n    \"\"\"crossover that selects random genes from the two parents\"\"\"\n    child = []\n    for p1, p2 in zip(parent1, parent2):\n        if np.random.random() < crossover_impact:\n            child.append(p1)\n        else:\n            child.append(p2)\n    return child\n\n\n# naive mutation\ndef mutate(individual):\n    \"\"\"mutate an individual by adding a random normal variable to each gene\"\"\"\n    for i in range(18):\n        if np.random.random() < mutation_impact:\n            individual[i] += np.random.normal(scale=0.001)\n    return individual\n\n\n# intelligent mutation\ndef mutate_intelligently(individual):\n    \"\"\"Mutation that shifts the thresholds towards separated thresholds\"\"\"\n    for i in range(18):\n        if np.random.random() < mutation_impact:\n            if individual[i] < thresholds_separated[i]:\n                individual[i] -= abs(np.random.normal(scale=0.001))\n            else:\n                individual[i] += abs(np.random.normal(scale=0.001))\n    return individual\n\n\n# Tournament selection\ndef tournament_selection(population, tournament_size, fitness_scores):\n    \"\"\"selection using tournament\"\"\"\n    tournament_indices = np.random.choice(len(population), size=tournament_size, replace=False)\n    tournament_fitness_scores = [fitness_scores[idx] for idx in tournament_indices]\n    max_fitness_index = np.argmax(tournament_fitness_scores)\n    return tournament_indices[max_fitness_index]\n\n\n# Genetic algorithm loop\nbest_score = 0\nno_improvement = 0\nfor generation in range(num_generations):\n    print('Generation', generation)\n    # Evaluate fitness for each individual\n    # Create a pool of worker processes with a maximum cores\n    with Pool(processes=os.cpu_count()) as pool:\n        # Evaluate fitness scores for each individual in the population using the process pool\n        fitness_scores = pool.map(evaluate_fitness, population)\n\n    # Find the best individual in the population\n    max_fitness_index = np.argmax(fitness_scores)\n    if fitness_scores[max_fitness_index] > best_score:\n        print('-------------------- best score improved --------------------')\n        best_score = fitness_scores[max_fitness_index]\n        best_thresholds = population[max_fitness_index]\n        no_improvement = 0\n    else:\n        if no_improvement > early_stopping:\n            print('-------------------- early stopping --------------------')\n            break\n        no_improvement += 1\n    best_scores.append(best_score)\n    print('best score:', best_score)\n    # calculate rmse of best score to original\n    print('rmse:', math.sqrt(mean_squared_error(best_thresholds, original_thresholds)))\n    print('best thresholds:', best_thresholds)\n    # write best thresholds to file\n    with open(output_file, 'a') as f:\n        f.write(f'{best_score},{math.sqrt(mean_squared_error(best_thresholds, original_thresholds))},{best_thresholds}\\n')\n\n    # Selection\n    selected_indices = []\n    for _ in range(population_size):\n        selected_indices.append(tournament_selection(population, tournament_size, fitness_scores))\n\n    # Crossover and mutation\n    offspring = []\n    for i in range(population_size):\n        parent1 = population[selected_indices[i]]\n\n        # Apply crossover with a probability\n        if np.random.random() < crossover_rate:\n            parent2 = population[tournament_selection(selected_indices, tournament_size, fitness_scores)]\n            if np.random.random() < 0.5:\n                child = average_crossover(parent1, parent2)\n            else:\n                child = selective_crossover(parent1, parent2)\n        else:\n            child = parent1\n\n        # Apply mutation with a probability\n        if np.random.random() < mutation_rate / 2:\n            child = mutate_intelligently(child)\n        if np.random.random() < mutation_rate / 2:\n            child = mutate(child)\n\n        offspring.append(child)\n\n    # Replace some individuals with the new offspring\n    elite_indices = np.argsort(fitness_scores)[-elite_size:]\n    for i in range(elite_size):\n        offspring[i] = population[elite_indices[i]]\n    population = offspring","metadata":{"execution":{"iopub.status.busy":"2023-06-02T11:41:49.074097Z","iopub.execute_input":"2023-06-02T11:41:49.074725Z","iopub.status.idle":"2023-06-02T12:46:46.840804Z","shell.execute_reply.started":"2023-06-02T11:41:49.074698Z","shell.execute_reply":"2023-06-02T12:46:46.839607Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# plot graph of improvement\nplt.figure(figsize=(12, 8))\nplt.title('Best score per generation')\nplt.xlabel('Generation')\nplt.ylabel('Score')\nplt.plot(best_scores)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-02T12:46:46.842314Z","iopub.execute_input":"2023-06-02T12:46:46.842619Z","iopub.status.idle":"2023-06-02T12:46:47.072106Z","shell.execute_reply.started":"2023-06-02T12:46:46.842590Z","shell.execute_reply":"2023-06-02T12:46:47.070736Z"},"trusted":true},"execution_count":null,"outputs":[]}]}