Genome_VAE.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "collapsed_sections": [
        "NAUaBLuBh3Fa"
      ]
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "SNQafZW7P4ut"
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "from IPython.display import display, HTML\n",
        "import pandas as pd\n",
        "from collections import Counter\n",
        "import plotly.express as px\n",
        "from tensorflow import keras\n",
        "from numba import jit\n",
        "from scipy.ndimage import gaussian_filter1d\n",
        "from scipy.signal import argrelmin\n",
        "from scipy.stats import entropy\n",
        "import tensorflow as tf\n",
        "from sklearn.model_selection import train_test_split\n",
        "import os\n",
        "from tensorflow.keras import layers\n",
        "from sklearn.metrics import confusion_matrix, accuracy_score\n",
        "import seaborn as sns\n",
        "import pickle\n",
        "from tqdm import tqdm\n",
        "import math\n",
        "from joblib import Parallel, delayed\n",
        "import multiprocessing\n",
        "from concurrent.futures import ThreadPoolExecutor\n",
        "from tensorflow.keras import regularizers\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.ndimage import gaussian_filter1d\n",
        "from scipy.signal import argrelmin\n",
        "\n",
        "\n",
        "print(\"Done!\")\n"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Load data\n",
        "---\n",
        "\n"
      ],
      "metadata": {
        "id": "ToPHiRxAQh3m"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#Example: Load data frame with Pickle\n",
        "path_to_data = \"C:/Users/ ..... /df_chr22_012.pkl\"\n",
        "\n",
        "with open(path_to_data, 'rb') as f:\n",
        "    df_012 = pickle.load(f)"
      ],
      "metadata": {
        "id": "VTEDLgKJQdVX"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#If no real data is available use this test data - You will however, get better results with real data\n",
        "num_rows = 1000\n",
        "num_cols = 100\n",
        "\n",
        "# Erstellen eines DataFrames mit zufälligen Nullen, Einsen und Zweien\n",
        "data = np.random.choice([0, 1, 2], size=(num_rows, num_cols))\n",
        "df_012 = pd.DataFrame(data, columns=[f'Individuum {i}' for i in range(num_cols)], index=[f'SNP {i}' for i in range(num_rows)])\n"
      ],
      "metadata": {
        "id": "DYqINm6IQ85R"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "\n",
        "Have a look at the data\n",
        "---\n",
        "\n"
      ],
      "metadata": {
        "id": "E06YHeZIR3lc"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "df_012"
      ],
      "metadata": {
        "id": "2ozcWolGR0mS"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Method 1:Split data frame into same size parts\n",
        "---\n",
        "\n"
      ],
      "metadata": {
        "id": "319P1HQSSDnG"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def split_dataframe(df, block_length):\n",
        "    num_blocks = len(df) // block_length\n",
        "    remainder = len(df) % block_length\n",
        "    block_indices = np.arange(0, num_blocks * block_length, block_length)\n",
        "\n",
        "\n",
        "    blocks = [df.iloc[i:i + block_length] for i in block_indices]\n",
        "\n",
        "    # if the are remains thy became there own segment\n",
        "    if remainder:\n",
        "        blocks.append(df.iloc[num_blocks * block_length:])\n",
        "\n",
        "    return blocks"
      ],
      "metadata": {
        "id": "LBMFts8ySD7J"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#Split data frame into parts of size 100\n",
        "segment_size = 100\n",
        "\n",
        "dataframes = split_dataframe(df_012,segment_size)\n",
        "total_length = 0\n",
        "for i in range(len(dataframes)):\n",
        "  print(\"Segment \",i,\" has a lenght of \",len(dataframes[i]))\n",
        "  total_length = total_length +len(dataframes[i])\n",
        "\n",
        "print(\"\")\n",
        "print(\"With a total length of\", total_length)"
      ],
      "metadata": {
        "id": "jnmhiOZYSapt"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Method 2: Split by correlation\n",
        "---\n",
        "\n"
      ],
      "metadata": {
        "id": "NAUaBLuBh3Fa"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#Generat correlation matrix\n",
        "@jit(parallel=True)\n",
        "def compute_correlation_matrix(data):\n",
        "    return np.corrcoef(data)\n",
        "\n",
        "\n",
        "correlation_matrix = compute_correlation_matrix(df_012.values)\n",
        "\n",
        "\n",
        "plt.figure(figsize=(10, 10))\n",
        "cmap = plt.cm.RdBu\n",
        "\n",
        "plt.imshow(correlation_matrix, cmap=cmap, vmin=-1, vmax=1)\n",
        "plt.colorbar()\n",
        "\n",
        "\n",
        "plt.title('Correlation matrix of variants from chromosome __', fontsize=20)\n",
        "plt.xlabel('Variant index', fontsize=16)\n",
        "plt.ylabel('Variant index', fontsize=16)\n",
        "\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "s48CkaNKhWCH"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#calculate aggregated correlation for segmentation\n",
        "\n",
        "threshould = 0.1 # this is the noise floor we found works fine for real data\n",
        "\n",
        "def antidiagonal_sum(matrix, row, col):\n",
        "    size = len(matrix)\n",
        "    diagonal_sum = 0\n",
        "\n",
        "    for i in range(size):\n",
        "        j = row + col - i\n",
        "        if j >= 0 and j < size:\n",
        "            if np.abs(matrix[i][j]) > threshould:\n",
        "                diagonal_sum += np.abs(matrix[i][j])\n",
        "\n",
        "    return diagonal_sum\n",
        "\n",
        "def calculate_antidiagonal_sums(matrix):\n",
        "    results = Parallel(n_jobs=-1)(delayed(antidiagonal_sum)(matrix, i, i) for i in range(len(matrix)))\n",
        "    return results\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "results = calculate_antidiagonal_sums(correlation_matrix)"
      ],
      "metadata": {
        "id": "klNbU5XNir0J"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#plot result\n",
        "\n",
        "#artificial data produces a pyramid , this nicely illustrates why we use a threshold for the real data with a threshold\n",
        "\n",
        "plt.figure(figsize=(16, 4))\n",
        "plt.plot(results)\n",
        "plt.xlabel('SNP Index')\n",
        "plt.ylabel('aggregated correlation')\n",
        "plt.title('raw data')\n",
        "plt.grid(True)\n",
        "plt.show()\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "smoothed_result = gaussian_filter1d(results, sigma=1)\n",
        "\n",
        "plt.figure(figsize=(16, 4))\n",
        "plt.plot(smoothed_result)\n",
        "\n",
        "plt.xlabel('SNP Index')\n",
        "plt.ylabel('aggregated correlation')\n",
        "plt.title('After gaussian filter')\n",
        "plt.grid(True)\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "MSPt785Akv6m"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#find local minimum\n",
        "\n",
        "\n",
        "\n",
        "smoothed_result = gaussian_filter1d(results, sigma=1)\n",
        "\n",
        "local_minima_indices = argrelmin(smoothed_result)[0]\n",
        "\n",
        "\n",
        "plt.figure(figsize=(16, 4))\n",
        "\n",
        "plt.plot(smoothed_result, label='Smoothed')\n",
        "plt.scatter(local_minima_indices, smoothed_result[local_minima_indices], color='red', label='Local Minima')\n",
        "plt.title('potential points of segmentation')\n",
        "plt.xlabel('SNP Index')\n",
        "plt.ylabel('aggregated correlation')\n",
        "plt.legend()\n",
        "plt.grid(True)\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "qSJJLMuMjqUb"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#Find best points for segmentation\n",
        "\n",
        "segment_size_min = 50#400\n",
        "segment_size_max = 150#1600\n",
        "\n",
        "position = 0\n",
        "\n",
        "segment_positions = []\n",
        "split_list = []\n",
        "\n",
        "run = True\n",
        "j = 0\n",
        "total_length_covered = 0\n",
        "while run:\n",
        "    candidates = []\n",
        "\n",
        "    for i in range(len(local_minima_indices)):\n",
        "        if local_minima_indices[i] >= (position + segment_size_min):\n",
        "            if local_minima_indices[i] <= (position + segment_size_max):\n",
        "                candidates.append(local_minima_indices[i])\n",
        "\n",
        "    if len(candidates) == 0:\n",
        "        break\n",
        "\n",
        "    minimum_x_position = candidates[np.argmin(smoothed_result[candidates])]\n",
        "    minimum_x_position_y_value = smoothed_result[minimum_x_position]\n",
        "\n",
        "    segment_positions.append([position, minimum_x_position])\n",
        "    print(\"Segment Nr:\", j)\n",
        "    print(\"Segment Start/End:\", [position, minimum_x_position])\n",
        "    print(\"Segment length:\", minimum_x_position - position)\n",
        "    total_length_covered += minimum_x_position - position\n",
        "    print(\"Total length covered:\", total_length_covered)\n",
        "    print(\"\")\n",
        "    split_list.append(minimum_x_position)\n",
        "    position = minimum_x_position + 1\n",
        "    j += 1\n",
        "\n",
        "    if position >= len(smoothed_result):\n",
        "        run = False\n",
        "\n",
        "\n",
        "if position < len(smoothed_result) and (len(smoothed_result) - position) < segment_size_min:\n",
        "    segment_positions.append([position, len(smoothed_result) - 1])\n",
        "    print(\"Segment Nr:\", j)\n",
        "    print(\"Segment Start/End:\", [position, len(smoothed_result) - 1])\n",
        "    print(\"Segment length:\", len(smoothed_result) - 1 - position)\n",
        "    total_length_covered += len(smoothed_result) - 1 - position\n",
        "    print(\"Total length covered:\", total_length_covered)"
      ],
      "metadata": {
        "id": "2r-1d8xPjz5b"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# plot the segmentation\n",
        "plt.figure(figsize=(30, 6))\n",
        "\n",
        "plt.plot(smoothed_result )\n",
        "\n",
        "for x_coord in split_list:\n",
        "    plt.axvline(x=x_coord, color='r', linestyle='--', linewidth=1 )\n",
        "\n",
        "\n",
        "plt.title('Segmentation plot')\n",
        "plt.xlabel('SNP index')\n",
        "plt.ylabel('aggregated correlation')\n",
        "#plt.legend()\n",
        "plt.grid(True)\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "4R-otxRzj5KG"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import pandas as pd\n",
        "\n",
        "def split_dataframe_by_positions(df, positions):\n",
        "    \"\"\"\n",
        "    Teilt ein DataFrame zeilenweise an den angegebenen Positionen auf.\n",
        "\n",
        "    Args:\n",
        "        df (pandas.DataFrame): Das DataFrame, das aufgeteilt werden soll.\n",
        "        positions (list): Eine Liste von Positionen, an denen das DataFrame aufgeteilt werden soll.\n",
        "\n",
        "    Returns:\n",
        "        list von pandas.DataFrame: Eine Liste von DataFrames, die durch das Aufteilen des ursprünglichen DataFrames an den angegebenen Positionen entstanden sind.\n",
        "    \"\"\"\n",
        "    # Sortiere die Positionen, um sicherzustellen, dass sie aufsteigend sind\n",
        "    positions.sort()\n",
        "\n",
        "    # Initialisiere die Liste für die aufgeteilten DataFrames\n",
        "    splitted_dfs = []\n",
        "\n",
        "    # Startindex für das Aufteilen des DataFrames\n",
        "    start_index = 0\n",
        "\n",
        "    # Iteriere über die Positionen und teile das DataFrame entsprechend auf\n",
        "    for pos in positions:\n",
        "        # Füge den aufgeteilten Teil des DataFrames zur Liste hinzu\n",
        "        splitted_dfs.append(df.iloc[start_index:pos])\n",
        "\n",
        "        # Setze den Startindex für den nächsten Teil\n",
        "        start_index = pos\n",
        "\n",
        "    # Füge den letzten Teil des DataFrames zur Liste hinzu\n",
        "    splitted_dfs.append(df.iloc[start_index:])\n",
        "\n",
        "    return splitted_dfs"
      ],
      "metadata": {
        "id": "x2NEr2g5j8CS"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#Split the data frame\n",
        "dataframes =  []\n",
        "dataframes = split_dataframe_by_positions(df_012,split_list)\n",
        "\n",
        "total_lenght_all_segments = 0\n",
        "start = 0\n",
        "\n",
        "for dataframe in range(len(dataframes)):\n",
        "\n",
        "  total_lenght_all_segments = total_lenght_all_segments +len(dataframes[dataframe])\n",
        "  #print(total_lenght_all_segments)\n",
        "  print(len(dataframes[dataframe]))\n",
        "\n",
        "print(\"total_lenght_all_segments:\",total_lenght_all_segments)"
      ],
      "metadata": {
        "id": "m8TzcDqmj_Y3"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Train VAE\n",
        "---\n",
        "\n"
      ],
      "metadata": {
        "id": "vKCmFx9QSuv3"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "\n",
        "def leaky_relu(x):\n",
        "    return tf.nn.leaky_relu(x, alpha=0.2)\n",
        "\n",
        "# Custom Aktivierungsfunktion r(x)\n",
        "def custom_activation(x):\n",
        "    return tf.math.tanh(x) + 1\n",
        "\n",
        "\n",
        "sum_all_variants = 0\n",
        "for elem in range(len(dataframes)):\n",
        "  sum_all_variants = sum_all_variants + dataframes[elem].shape[0]\n",
        "\n",
        "\n",
        "reconstruction = 0\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "for Block_ID in range(len(dataframes)):\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "  data = dataframes[Block_ID].T.values\n",
        "\n",
        "# Hyperparameter\n",
        "######################################################################################################################################################################################\n",
        "######################################################################################################################################################################################\n",
        "######################################################################################################################################################################################\n",
        "  train_val_ratio = 0.8#  the ratio of how large the training set should be\n",
        "\n",
        "\n",
        "  # splitting into test and validation data\n",
        "  train_size = int(train_val_ratio * len(data))\n",
        "  train_data = data[:train_size]\n",
        "  val_data = data[train_size:]\n",
        "\n",
        "\n",
        "  foldername = \"Run\"                                      #in case one wanst to have multiple training runs you can add the iteration to this name:   foldername = \"Run\"   + str(iteration_nr)\n",
        "  experiment_name = \"VAE_Demo\"                            #name of the experiment\n",
        "  bottleneck_size = max(1, int(train_data.shape[1]*0.1))   #how large the embedding should be in relationship to the input data:  currently 10%\n",
        "  latent_dim = bottleneck_size\n",
        "  batch_size = 256\n",
        "  epochs = 5000                                           #the maximum number of epochs for training\n",
        "  patience = 30                                           #after this many epox of no improvement the training will be terminated\n",
        "  chromosome = 22                                       #part of the folder name structure\n",
        "  dropout = 0.5                                           #drop outvalue of the dropout layers\n",
        "  activation_func = leaky_relu\n",
        "  learn_rate = 1e-4\n",
        "  gradient_clip_value = 0.25\n",
        "  fixed_beta = 0.015                                      #currently using fixed beta value\n",
        "  output_steps = 0                                        # output additional information during training:  0 off , 1 on\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "  alpha = 1\n",
        "  beta = 1\n",
        "\n",
        "  Start_KLAnnealing_Epoch = 1\n",
        "  Start_value_KLAnnealing = fixed_beta\n",
        "  epoch_history = []\n",
        "  beta_history = []\n",
        "\n",
        "  #this works for online Google collab but change if you want to use locally\n",
        "  stotage_path =  \"C:/Users/HD-ThinkTank/Desktop/Chr\" + str(chromosome) + \"_models/\" + experiment_name + \"/\"+foldername+ \"/\" + \"Block_ID_\" + str(Block_ID) +\"/\"\n",
        "  stotage_path_txt =    \"C:/Users/HD-ThinkTank/Desktop/Chr\" + str(chromosome) +  \"_models/\" + experiment_name +  \"/\"+foldername+ \"/\"\n",
        "\n",
        "#####No_Settings_after_this_point##########No_Settings_after_this_point##########No_Settings_after_this_point##########No_Settings_after_this_point##########No_Settings_after_this_point#####\n",
        "######################################################################################################################################################################################\n",
        "######################################################################################################################################################################################\n",
        "######################################################################################################################################################################################\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "  print(\"##############################################\")\n",
        "  print(\"Now training block \", Block_ID,\"/\",len(dataframes)-1)\n",
        "  print(\"bottleneck_size: \",bottleneck_size)\n",
        "  print(\"Variants: \", train_data.shape[1] )\n",
        "  print(\"datapoints: \",train_data.shape[0])\n",
        "  print(\"compression:  \",train_data.shape[1]/bottleneck_size)\n",
        "  print(\"##############################################\")\n",
        "\n",
        "  # Variational Autoencoder Klasse\n",
        "  class VAE(keras.Model):\n",
        "      def __init__(self, latent_dim, dropout_rate=dropout):\n",
        "          super(VAE, self).__init__()\n",
        "          self.latent_dim = latent_dim\n",
        "          self.dropout_rate = dropout_rate\n",
        "          self.beta = tf.Variable(0.0, trainable=False, dtype=tf.float32)  # Initial beta value\n",
        "\n",
        "\n",
        "          self.encoder = tf.keras.Sequential([\n",
        "              #Input Layer\n",
        "              layers.InputLayer(input_shape=(data.shape[1],)),\n",
        "\n",
        "              #Hidden Block\n",
        "              layers.Dense(int(data.shape[1]), activation=activation_func),\n",
        "              layers.BatchNormalization(),  # Batch-Normalization-\n",
        "              layers.Dropout(dropout_rate),\n",
        "\n",
        "              #Hidden Block\n",
        "              layers.Dense(int(data.shape[1]*0.5), activation=activation_func),\n",
        "              layers.BatchNormalization(),  # Batch-Normalization-\n",
        "              layers.Dropout(dropout_rate),\n",
        "\n",
        "              #Output Layer\n",
        "              layers.Dense(latent_dim+latent_dim)\n",
        "          ])\n",
        "\n",
        "          self.decoder = tf.keras.Sequential([\n",
        "              #Input Layer\n",
        "              layers.InputLayer(input_shape=(latent_dim,)),\n",
        "\n",
        "              #Hidden Block\n",
        "              layers.Dense(int(data.shape[1]*0.5), activation=activation_func),\n",
        "              layers.BatchNormalization(),  # Batch-Normalization\n",
        "\n",
        "              #Hidden Block\n",
        "              layers.Dense(int(data.shape[1]), activation=activation_func),\n",
        "              layers.BatchNormalization(),  # Batch-Normalization\n",
        "\n",
        "              #Output Layer\n",
        "              layers.Dense(data.shape[1], activation=custom_activation)\n",
        "          ])\n",
        "\n",
        "      def encode(self, x):\n",
        "          mean, logvar = tf.split(self.encoder(x), num_or_size_splits=2, axis=1)\n",
        "          return mean, logvar\n",
        "\n",
        "      def reparameterize(self, mean, logvar):\n",
        "          eps = tf.random.normal(shape=tf.shape(mean))\n",
        "          return eps * tf.exp(logvar * .5) + mean\n",
        "\n",
        "      def decode(self, z):\n",
        "          return self.decoder(z)\n",
        "\n",
        "      def call(self, x):\n",
        "          mean, logvar = self.encode(x)\n",
        "          z = self.reparameterize(mean, logvar)\n",
        "          x_recon = self.decode(z)\n",
        "          return x_recon, mean, logvar\n",
        "\n",
        "  # VAE Modell erstellen\n",
        "  vae = VAE(latent_dim, dropout)\n",
        "\n",
        "  # Optimizer definieren mit Gradient Clipping und einstellbarer Lernrate\n",
        "  optimizer = tf.keras.optimizers.Adam(learning_rate=learn_rate, clipvalue=gradient_clip_value)  # Hier kannst du den Clip-Wert anpassen\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "  # KL-Annealing-Callback definieren\n",
        "  class KLAnnealingCallback(tf.keras.callbacks.Callback):\n",
        "    def __init__(self, beta_start= fixed_beta, beta_end=fixed_beta, n_epochs=Start_KLAnnealing_Epoch):\n",
        "        super(KLAnnealingCallback, self).__init__()\n",
        "        self.beta_start = beta_start\n",
        "        self.beta_end = beta_end\n",
        "        self.n_epochs = n_epochs\n",
        "\n",
        "    def on_epoch_begin(self, epoch, logs=None):\n",
        "        if epoch < self.n_epochs:\n",
        "            new_beta = fixed_beta\n",
        "            vae.beta.assign(new_beta)\n",
        "        else:\n",
        "            #vae.beta.assign(self.beta_end)\n",
        "            vae.beta.assign(fixed_beta)\n",
        "            new_beta = fixed_beta# self.beta_end\n",
        "\n",
        "        beta_history.append(new_beta)\n",
        "        epoch_history.append(epoch)\n",
        "\n",
        "  def vae_loss(x, x_recon):\n",
        "    x_recon, mean, logvar = vae(x)\n",
        "    recon_loss = tf.keras.losses.MeanSquaredError()(x, x_recon)\n",
        "    kl_loss = -0.5 * tf.reduce_mean(1 + logvar - tf.square(mean) - tf.exp(logvar), axis=-1)\n",
        "    x = tf.cast(x, tf.float32)\n",
        "    #tf.print(\"with a beta of \", vae.beta )\n",
        "    return  recon_loss + vae.beta * kl_loss + 1e-8\n",
        "\n",
        "  # Modell Checkpoint für die besten Gewichte speichern\n",
        "  checkpoint_filepath = stotage_path\n",
        "  model_checkpoint_callback = tf.keras.callbacks.ModelCheckpoint(\n",
        "      filepath=checkpoint_filepath,\n",
        "      save_weights_only=True,\n",
        "      save_best_only=True,\n",
        "      monitor='val_loss',\n",
        "      mode='min',\n",
        "      verbose=0\n",
        "  )\n",
        "\n",
        "  # Modell kompilieren\n",
        "  vae.compile(optimizer=optimizer, loss=vae_loss)\n",
        "\n",
        "  class CustomEarlyStopping(tf.keras.callbacks.Callback):\n",
        "      def __init__(self, monitor='val_loss', patience=30, start_epoch=Start_KLAnnealing_Epoch):\n",
        "          super(CustomEarlyStopping, self).__init__()\n",
        "          self.monitor = monitor\n",
        "          self.patience = patience\n",
        "          self.start_epoch = start_epoch\n",
        "          self.wait = 0\n",
        "          self.best = float('inf')\n",
        "\n",
        "      def on_epoch_end(self, epoch, logs=None):\n",
        "          current = logs.get(self.monitor)\n",
        "          if epoch >= self.start_epoch:\n",
        "              if current < self.best:\n",
        "                  self.best = current\n",
        "                  self.wait = 0\n",
        "              else:\n",
        "                  self.wait += 1\n",
        "                  if self.wait >= self.patience:\n",
        "                      self.model.stop_training = True\n",
        "\n",
        "\n",
        "  # Modell trainieren mit Early Stopping Callback\n",
        "  history = vae.fit(train_data, train_data,\n",
        "                  batch_size=batch_size,\n",
        "                  epochs=epochs,\n",
        "                  validation_data=(val_data, val_data),\n",
        "                  callbacks=[model_checkpoint_callback, KLAnnealingCallback(beta_start=Start_value_KLAnnealing), CustomEarlyStopping()],\n",
        "                  verbose=output_steps)\n",
        "  #---------------------------------------------------------------------------------------------------\n",
        "  # Ausgabe des Verlaufs\n",
        "\n",
        "\n",
        "\n",
        "  plt.plot(epoch_history,beta_history, label='β value')\n",
        "  plt.plot(history.history['loss'], label='Training Loss')\n",
        "  plt.plot(history.history['val_loss'], label='Validation Loss')\n",
        "  plt.xlabel('Epoch')\n",
        "  plt.ylabel('Loss')\n",
        "  plt.legend()\n",
        "  plt.show()\n",
        "\n",
        "  # Bestes Modell laden\n",
        "  vae.load_weights(checkpoint_filepath)\n",
        "\n",
        "  #---------------------------------------------------------------------------------------------------\n",
        "\n",
        "  num_runs = 100\n",
        "  total_percentage_correct_reconstruction = 0\n",
        "\n",
        "  for _ in range(num_runs):\n",
        "      val_reconstructions, _, _ = vae(val_data)\n",
        "      mse = tf.keras.losses.MSE(val_data, tf.round(val_reconstructions))\n",
        "      mean_mse = tf.reduce_mean(mse)\n",
        "      percentage_correct_reconstruction = (1 - mean_mse) * 100\n",
        "      total_percentage_correct_reconstruction += percentage_correct_reconstruction\n",
        "\n",
        "  average_percentage_correct_reconstruction = total_percentage_correct_reconstruction / num_runs\n",
        "  print(f\"Avarage reconstruction accuracy over {num_runs} runs: {average_percentage_correct_reconstruction:.2f}%\")\n",
        "\n",
        "  with open(stotage_path +'Hyperparameter.txt', 'w') as file:\n",
        "      file.write(f\"bottleneck_size: {bottleneck_size}\\n\")\n",
        "      file.write(f\"Variablen: {train_data.shape[1]}\\n\")\n",
        "      file.write(f\"epochs: {epochs}\\n\")\n",
        "      file.write(f\"patience: {patience}\\n\")\n",
        "      file.write(f\"Block_ID: {Block_ID}\\n\")\n",
        "      file.write(f\"average_percentage_correct_reconstruction: {average_percentage_correct_reconstruction}\\n\")\n",
        "\n",
        "  reconstruction = reconstruction + float(average_percentage_correct_reconstruction * (train_data.shape[1]/sum_all_variants ))\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "  # Funktion zur Berechnung der KL-Divergenz\n",
        "  def compute_kl_divergence(mean, logvar):\n",
        "      return -0.5 * np.mean(1 + logvar - np.square(mean) - np.exp(logvar), axis=-1)\n",
        "\n",
        "  # Berechnen der KL-Divergenz für die Validierungsdaten\n",
        "  def evaluate_kl_divergence(model, val_data):\n",
        "      mean, logvar = model.encode(val_data)\n",
        "      kl_divergence = compute_kl_divergence(mean.numpy(), logvar.numpy())\n",
        "      return np.mean(kl_divergence)\n",
        "\n",
        "  # Berechnen und ausgeben der KL-Divergenz für die Validierungsdaten\n",
        "  kl_divergence = evaluate_kl_divergence(vae, val_data)\n",
        "  print(f\"KL Divergence on validation data: {kl_divergence}\")\n",
        "\n",
        "  with open(stotage_path +'kl_divergence.txt', 'w') as file:\n",
        "    file.write(f\"reconstruction: {kl_divergence}\")\n",
        "\n",
        "\n",
        "\n",
        "  ####-----------------------------------\n",
        "\n",
        "\n",
        "  # Funktion zur Berechnung der KL-Divergenz\n",
        "  def compute_kl_divergence(mean, logvar):\n",
        "      return -0.5 * (1 + logvar - np.square(mean) - np.exp(logvar))\n",
        "\n",
        "  # Berechnen der KL-Divergenz für die Validierungsdaten\n",
        "  def evaluate_kl_divergence(model, val_data):\n",
        "      mean, logvar = model.encode(val_data)\n",
        "      kl_divergence_per_dim = compute_kl_divergence(mean.numpy(), logvar.numpy())\n",
        "      kl_divergence_mean = np.mean(kl_divergence_per_dim, axis=0)\n",
        "      kl_divergence_total = np.sum(kl_divergence_mean)\n",
        "      return kl_divergence_total, kl_divergence_mean\n",
        "\n",
        "  # Berechnen und ausgeben der KL-Divergenz für die Validierungsdaten\n",
        "  kl_divergence_total, kl_divergence_per_dim = evaluate_kl_divergence(vae, val_data)\n",
        "  print(f\"KL Divergence on validation data (total): {kl_divergence_total}\")\n",
        "  print(f\"KL Divergence on validation data (per dimension): {kl_divergence_per_dim}\")\n",
        "\n",
        "  # Berechnen und ausgeben der KL-Divergenz für die Validierungsdaten\n",
        "  kl_divergence = evaluate_kl_divergence(vae, val_data)\n",
        "  print(f\"KL Divergence on validation data: {kl_divergence}\")\n",
        "\n",
        "  # Daten durch den Encoder leiten, um die latenten Mittelwerte und Log-Varianzen zu erhalten\n",
        "  mean, logvar = vae.encode(val_data)\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "def ordner_erstellen(pfad):\n",
        "  try:\n",
        "      os.makedirs(pfad)\n",
        "      print(\"Ordner erfolgreich erstellt:\", pfad)\n",
        "  except FileExistsError:\n",
        "      print(\"Der Ordner existiert bereits.\")\n",
        "  except Exception as e:\n",
        "      print(\"Fehler beim Erstellen des Ordners:\", e)\n",
        "\n",
        "# Beispielaufruf der Funktion\n",
        "ordner_erstellen(stotage_path_txt)\n",
        "\n",
        "with open(stotage_path_txt +'reconstruction.txt', 'w') as file:\n",
        "    file.write(f\"reconstruction: {reconstruction}\")\n",
        "\n"
      ],
      "metadata": {
        "id": "1oCh-Bm0Schd"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}