Similarity_QUBO_processing.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pickle\n",
    "import dimod\n",
    "import os\n",
    "import dwave_token\n",
    "import dwave.inspector\n",
    "\n",
    "from dwave.system import DWaveSampler, EmbeddingComposite, FixedEmbeddingComposite\n",
    "from dwave.system import DWaveCliqueSampler\n",
    "# from dwave.system import LeapHybridCQMSampler\n",
    "\n",
    "# automatically generated embedding:\n",
    "sampler = EmbeddingComposite(DWaveSampler(token=dwave_token.value))\n",
    "clique_sampler = DWaveCliqueSampler(token=dwave_token.value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def csv_to_matrix(file_name): # read csv file to get row names and matrix values\n",
    "    dimension = len(np.genfromtxt(file_name)) - 1\n",
    "    rows = np.array([element[0] for element in pd.read_csv(file_name, delimiter=';', usecols=[0]).to_numpy()])\n",
    "    matrix = np.genfromtxt(file_name, delimiter=';', skip_header=1, usecols = range(1,dimension+1))\n",
    "    return(rows, matrix)\n",
    "\n",
    "def distance_to_similarity(Q, neutral_distance=1): # convert distance matrix to similarity\n",
    "    # neutral distance is the distance which will be mapped to similarity 0, larger distance values will be mapped to negative similarity\n",
    "    # The worst similarity (distance 1) will then be (neutral_distance-1)/neutral_distance\n",
    "    P = 1 - np.array(Q)/neutral_distance\n",
    "    return(P)\n",
    "\n",
    "def load_mapping(file_name): # load mapping file and convert to normalized number values\n",
    "    frame = pd.read_csv(file_name, delimiter=';')\n",
    "    class_to_number = {\n",
    "        \"Approved\": 1,\n",
    "        \"Preclinical\": 0,\n",
    "        \"Withdrawn\": -1\n",
    "    }\n",
    "    dictionary = {key: class_to_number[val] for key, val in frame.set_index('ROWID')['Class'].to_dict().items()}\n",
    "    return(frame, dictionary)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\max \\left(\\sum_i h_i \\cdot s_i + \\sum_{i<j} J_{i,j} \\cdot s_i \\cdot s_j \\right)$\n",
    "\n",
    "Linear bias of compounds:\n",
    "$$h_i = \n",
    "    \\begin{cases}\n",
    "        \\text{bias\\_scaling} &, \\text{compound $i$ approved}\\\\\n",
    "        0 &, \\text{compound $i$ preclinical}\\\\\n",
    "        - \\text{bias\\_scaling} &, \\text{compound $i$ failed/withdrawn}\n",
    "    \\end{cases}\n",
    "$$ \n",
    "Quadratic bias:\n",
    "$$J_{i,j} = 1 - \\text{distance}_{i,j}/\\text{neutral\\_distance}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "folder = 'Compound data/'\n",
    "file_name = 'Matrix_ECFP4.csv'\n",
    "mapping_file_name = 'mapping.csv'\n",
    "N = len(csv_to_matrix(os.path.join(folder,file_name))[0])\n",
    "bias_scaling = 0.05*N # choose scaling for linear bias\n",
    "neutral_distance = 0.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "rows, distance_J = csv_to_matrix(os.path.join(folder,file_name)) # row names and distances\n",
    "mapping, bias_dict0 = load_mapping(os.path.join(folder, mapping_file_name)) # linear biases from the mapping file\n",
    "J = np.triu(distance_to_similarity(distance_J,neutral_distance=neutral_distance)) # convert into a similarity upper triangle matrix. set distance that is without reward/penalty\n",
    "J_dict = {(rows[i],rows[j]): -J[i,j] for i in range(N) for j in range(i+1,N)} # - for minimization\n",
    "bias = np.array([bias_scaling * bias_dict0[row] for row in rows])\n",
    "bias_dict = {key: -bias_scaling * value for (key,value) in bias_dict0.items()} # - for minimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # if feasible, load embedding:\n",
    "# embedding_folder = 'Embeddings' # folder where embedding is saved\n",
    "# #embedding_folder = folder\n",
    "# embedding_file_name = str(N) + '_Ising' + '_embedding.pkl'\n",
    "# #embedding_file_name =  os.path.splitext(file_name)[0] + '_embedding.pkl'\n",
    "# with open(os.path.join(embedding_folder,embedding_file_name), 'rb') as f:\n",
    "#     embedding = pickle.load(f)\n",
    "# if len(embedding) != N:\n",
    "#     print('Choose and embedding of correct size!')\n",
    "# elif rows[0] not in embedding.keys():\n",
    "#     try:\n",
    "#         embedding = {row: embedding[i] for i, row in enumerate(rows)}\n",
    "#     except KeyError:\n",
    "#         embedding = {row: list(embedding.values())[i] for i, row in enumerate(rows)}\n",
    "# fixed_sampler = FixedEmbeddingComposite(child_sampler=DWaveSampler(token=dwave_token.value), embedding=embedding)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "bqm = dimod.BinaryQuadraticModel.from_ising(bias_dict,J_dict) # model creation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "minimum: -206.6341052883971\n",
      "102 variables encoded in 1002 physical qubits\n"
     ]
    }
   ],
   "source": [
    "chain_strength = bias_scaling # set chain strength close to largest coefficient\n",
    "# sampleset = sampler.sample( # uncomment this if you want to generate an embedding on the fly\n",
    "# sampleset = fixed_sampler.sample( # uncomment this if you want to used the fixed saved embedding\n",
    "sampleset = clique_sampler.sample( # uncomment this if you want to use a presaved embedding of dense problems\n",
    "    bqm,\n",
    "    num_reads=1000,\n",
    "    chain_strength=chain_strength\n",
    ")\n",
    "print('minimum: ' + str(sampleset.lowest().first.energy))\n",
    "print(str(N) + ' variables encoded in ' + str(len([node for nodes in sampleset.info['embedding_context']['embedding'].values() for node in nodes])) + ' physical qubits')\n",
    "solution_dict = sampleset.lowest().first.sample # load best solution in a dictionary\n",
    "solution = np.array([solution_dict[row] for row in rows])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Approved which are bad: 0\n",
      " Withdrawn which are good: 0\n",
      " Preclinical which are bad 44\n",
      " Preclinical which are good 7\n",
      "\n",
      "[ 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1\n",
      "  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1\n",
      "  1  1  1  1  1  1  1  1  1  1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1\n",
      " -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1\n",
      " -1 -1 -1 -1 -1 -1]\n"
     ]
    }
   ],
   "source": [
    "print(\n",
    "    'Approved which are bad: ' + str(np.count_nonzero(np.logical_and(np.sign(bias) == 1, solution < 0))) + '\\n',\n",
    "    'Withdrawn which are good: ' + str(np.count_nonzero(np.logical_and(np.sign(bias) == -1, solution > 0))) + '\\n',\n",
    "    'Preclinical which are bad ' + str(np.count_nonzero(np.logical_and(np.sign(bias) == 0, solution < 0))) + '\\n',\n",
    "    'Preclinical which are good ' + str(np.count_nonzero(np.logical_and(np.sign(bias) == 0, solution > 0))) + '\\n'\n",
    "    )\n",
    "print(solution)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save embedding:\n",
    "embedding_file_name = os.path.splitext(file_name)[0] + '_embedding.pkl'\n",
    "with open(os.path.join(folder,embedding_file_name), 'wb') as f:\n",
    "    pickle.dump(sampleset.to_serializable()['info']['embedding_context']['embedding'], f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparison with classical computation:\n",
    "Calculating similarity with approved compounds and then ordering accordingly\n",
    "\n",
    "$sim_i = \\sum_{j\\neq i\\ \\text{approved}} \\text{similarity}(i,j)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [],
   "source": [
    "similarity_dict = {(rows[i],rows[j]): distance_to_similarity(distance_J,0.3)[i,j] for i in range(N) for j in range(N)}\n",
    "approved_list = [row for row in rows if np.sign(bias_dict0[row]) == 1]\n",
    "sim_dict = {row: sum(similarity_dict[row,pow] for pow in approved_list if pow != row) for row in rows}\n",
    "sorted_dict = sorted(sim_dict.items(), key=lambda x:x[1], reverse=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[('Row25', -94.07054901472613), ('Row50', -94.794893548157), ('Row39', -96.35794760686892), ('Row32', -96.3900904668777), ('Row37', -96.69211663090114), ('Row22', -96.91378022072725), ('Row11', -96.9213318987318), ('Row0', -97.28615864743844), ('Row7', -97.5112293736308), ('Row2', -97.57237662470456)]\n"
     ]
    }
   ],
   "source": [
    "print(sorted_dict[:10])\n",
    "# It seems like with the neutral distance of 0.3, all compounds are far away from the group"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n",
      "True\n",
      "True\n",
      "True\n",
      "True\n",
      "True\n",
      "True\n",
      "True\n",
      "True\n",
      "True\n",
      "True\n",
      "True\n",
      "True\n",
      "True\n",
      "True\n",
      "True\n",
      "True\n",
      "False\n",
      "True\n",
      "True\n",
      "True\n",
      "False\n",
      "True\n",
      "True\n",
      "True\n",
      "False\n",
      "True\n",
      "True\n",
      "True\n",
      "True\n",
      "False\n",
      "False\n",
      "True\n",
      "True\n",
      "True\n",
      "False\n",
      "False\n",
      "False\n",
      "True\n",
      "False\n",
      "True\n",
      "True\n",
      "False\n",
      "False\n",
      "True\n",
      "False\n",
      "False\n",
      "False\n",
      "True\n",
      "False\n",
      "False\n",
      "True\n",
      "False\n",
      "True\n",
      "True\n",
      "False\n",
      "True\n",
      "False\n",
      "False\n",
      "False\n",
      "False\n",
      "False\n",
      "True\n",
      "False\n",
      "False\n",
      "False\n",
      "False\n",
      "False\n",
      "False\n",
      "False\n",
      "False\n",
      "True\n",
      "False\n",
      "True\n",
      "False\n",
      "False\n",
      "False\n",
      "False\n",
      "True\n",
      "True\n",
      "True\n",
      "False\n",
      "False\n",
      "False\n",
      "False\n",
      "True\n",
      "False\n",
      "False\n",
      "False\n",
      "False\n",
      "True\n",
      "True\n",
      "False\n",
      "False\n",
      "False\n",
      "False\n",
      "False\n",
      "False\n",
      "True\n",
      "True\n",
      "False\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "for key,val in sorted_dict:\n",
    "    print(key in approved_list)\n",
    "# But at least it looks like the approved group is the closest to itself"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}