samples/notebooks/ProteinFolding.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "dac00205-4c02-4e30-8bb2-f08b5ae7a200",
   "metadata": {},
   "source": [
    "# Simple Protein Folding Optimization\n",
    "\n",
    "This notebook is derived from the [Protein Folding](https://www.gurobi.com/jupyter_models/protein-folding/)\n",
    "Jupyter notebook example published by Gurobi.\n",
    "\n",
    "Here, the setting is the same, but we use Rexl rather than Python to formulate the model.\n",
    "Rexl includes support for multiple MIP solvers, including [Gurobi](https://www.gurobi.com/),\n",
    "[HiGHS](https://highs.dev/), and [GLPK](https://www.gnu.org/software/glpk/).\n",
    "This makes it easy to compare their relative performance."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "432ffad4-8e7c-4e5d-b361-ff246f4f64f6",
   "metadata": {},
   "source": [
    "## Problem Description\n",
    "\n",
    "For background, see the\n",
    "[source notebook](https://colab.research.google.com/github/Gurobi/modeling-examples/blob/master/protein_folding/protein_folding.ipynb)\n",
    "published by Gurobi.\n",
    "\n",
    "We summarize here. A certain protein consists of `Num` amino acids arranged linearly.\n",
    "A subset of the amino acids are _hydro-phobic_. The protein tends to fold so that\n",
    "the number of _matches_ of _hydro-phobic_ amino acids is maximized.\n",
    "\n",
    "A fold at position `R` allows positions `P` and `Q` to match if `R` and `R+1` are half\n",
    "way between `P` and `Q`. That is, a fold at `R` allows `P` and `Q` to match if `R = (P + Q - 1) div 2`.\n",
    "However, if there are any other folds between `P` and `Q` then they will not be matched.\n",
    "\n",
    "These pictures are taken from The Gurobi notebook. They show a protein in unfolded and folded forms,\n",
    "with the hydro-phobic amino acids shaded.\n",
    "\n",
    "![chain](https://github.com/Gurobi/modeling-examples/blob/master/protein_folding/chain.PNG?raw=1)\n",
    "\n",
    "![folding](https://github.com/Gurobi/modeling-examples/blob/master/protein_folding/folding.PNG?raw=1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9994b272-75fa-4389-9904-a269e7123d06",
   "metadata": {},
   "source": [
    "## Model Formulation\n",
    "\n",
    "Rexl supports the concept of module, which contains various kinds of symbols:\n",
    "* A _parameter_ (`param`) is defined by a formula. The formula can reference other parameters and constants.\n",
    "* A _constant_ (`const`) is defined by a formula. The formula can reference parameters and other constants.\n",
    "* A _free variable_ (`var`) is defined by its domain. There are multiple ways to specify a domain.\n",
    "  The formulas that define the domain can reference only parameters and constants, not variables.\n",
    "* A _computed variable_ is defined by a formula. The formula can reference any kind of symbol.\n",
    "  There are multiple kinds of computed variables:\n",
    "  * A measure (`msr`) is typically used as an objective (to maximize or minimize).\n",
    "  * A constraint (`con`) is boolean valued. A solution will make these be `true`.\n",
    "  * A let (`let`) is typically for displaying information or as an intermediate result used by other computed variables.\n",
    "\n",
    "New instances of a module can be constructed with different values for the parameters and free variables.\n",
    "When a solver is invoked to maximize or minimize a measure, the result (when successful) is a new instance\n",
    "of the module with the free variable values set."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e7b3821-126e-41a0-9df1-cc9ef8300681",
   "metadata": {},
   "source": [
    "We now define the `Folding` module. It has parameters for the number of positions in the protein\n",
    "and the hydro-phobic positions. Below, we encode an alternate form for the parameters, when the\n",
    "protein is specified via a text value containing `H` for a hydro-phobic location and `-` for others."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8452bc27-3935-43e4-bd05-a8ecc5a4cabc",
   "metadata": {},
   "outputs": [],
   "source": [
    "Folding := module {\n",
    "    // The number of amino acids in the protein.\n",
    "    // The positions of the hydro-phobic amino acids.\n",
    "    param Num := 50;\n",
    "    param HydroPhobic := [2,4,5,6,11,12,17,20,21,25,27,28,30,31,33,37,44,46];\n",
    "    const NumHydroPhobic := HydroPhobic->Count();\n",
    "\n",
    "    // Construct a table of possible matches, including the matching positions (P and Q)\n",
    "    // and the fold location (R). Then augment with an assigned index (Id) and with\n",
    "    // the sequence of bad fold locations (X).\n",
    "    const PossibleMatches :=\n",
    "        CrossJoin(P: HydroPhobic, Q: HydroPhobic,\n",
    "            Q - P >= 3 and (Q - P) mod 2 = 1, // The join predicate.\n",
    "            { P, Q, R: (P + Q - 1) div 2 })   // The result.\n",
    "        // Add the index (Id) and bad fold locations (X).\n",
    "        +>{ Id: #, X: Range(P, Q)->TakeIf(it != R) };\n",
    "\n",
    "    // The number of possible matches.\n",
    "    const NumPossibleMatches := PossibleMatches->Count();\n",
    "\n",
    "    // Define the free variables:\n",
    "    // * A true/false variable for each possible match.\n",
    "    // * A true/false variable for each possible fold location.\n",
    "    var Matches := Tensor.Fill(false, NumPossibleMatches);\n",
    "    var Folds := Tensor.Fill(false, Num);\n",
    "\n",
    "    // We want to maximize the number of matches.\n",
    "    msr NumMatches := Sum(Matches.Values);\n",
    "\n",
    "    // Construct the constraints.\n",
    "    // For a match, there must be a fold at R.\n",
    "    // For a match, there must NOT be a fold at any of the bad fold locations (X).\n",
    "    con MustFold := PossibleMatches->ForEach(Matches[Id] <= Folds[R]);\n",
    "    con CantFold := PossibleMatches->ForEach(Matches[Id] + Folds[X] <= 1);\n",
    "\n",
    "    // Define the table of matches.\n",
    "    let MatchTable := PossibleMatches->TakeIf(Matches[Id])->{ Id, P, Q, R };\n",
    "};"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "428c444b-5fed-42a0-8a0a-71ace18fe48f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "50\r\n",
       "18\r\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Folding.Num;\n",
    "Folding.NumHydroPhobic;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "001018cc-523a-4d65-b0a1-7d81dc8a59ab",
   "metadata": {},
   "source": [
    "## Optimize\n",
    "\n",
    "Now optimize using each of the supported MIP solvers.  We run each twice to see both cold and warm\n",
    "solve times. [HiGHS](https://highs.dev/) really shines on this problem.\n",
    "\n",
    "If you don't have a Gurobi license, the second one will fail (result in `null`). Also note that the first\n",
    "run of Gurobi includes time to verify the license, which may artificially inflate the time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e4771ffe-23d2-4395-bc5f-751215891fc4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Solver: HiGHS\r\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Wall time: 135.6336ms"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#!time\n",
    "Best1 := Folding->Maximize(NumMatches, \"highs\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c4f0bdff-aa4e-4099-a079-b8eaa0bc3801",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Solver: HiGHS\r\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Wall time: 72.9459ms"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#!time\n",
    "Best1 := Folding->Maximize(NumMatches, \"highs\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "84fcbd98-223f-45f4-ae02-051d826a0279",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Solver: Gurobi\r\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Wall time: 3039.4419ms"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#!time\n",
    "Best2 := Folding->Maximize(NumMatches, \"gurobi\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "7e12c2c8-bb9b-4c03-99c5-d6b2551ff237",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Solver: Gurobi\r\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Wall time: 327.4464ms"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#!time\n",
    "Best2 := Folding->Maximize(NumMatches, \"gurobi\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "5f680d3d-f3ac-4c1d-b4a2-e36f3fc30562",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Solver: GLPK\r\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Wall time: 2702.3877ms"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#!time\n",
    "Best3 := Folding->Maximize(NumMatches, \"glpk\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "d2109ca5-2999-4686-ace8-769f9c56650d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Solver: GLPK\r\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Wall time: 2687.5423ms"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#!time\n",
    "Best3 := Folding->Maximize(NumMatches, \"glpk\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "ee2ef144-d0a1-434c-883a-a134d5e38a10",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10\r\n",
       "10\r\n",
       "10\r\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Best1.NumMatches;\n",
    "Best2.NumMatches;\n",
    "Best3.NumMatches;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e27df7cf-f5fd-4e7a-83fe-8a1581a7b0e6",
   "metadata": {},
   "source": [
    "See if the solutions match."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "d73fbdad-bd90-4930-998c-c95a2acb4be0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Seq<{i8,i8,i8,i8}>\r\n",
       "   0) { Id: 0, P: 2, Q: 5, R: 3 }\r\n",
       "   1) { Id: 17, P: 5, Q: 12, R: 8 }\r\n",
       "   2) { Id: 23, P: 6, Q: 11, R: 8 }\r\n",
       "   3) { Id: 36, P: 12, Q: 17, R: 14 }\r\n",
       "   4) { Id: 43, P: 17, Q: 20, R: 18 }\r\n",
       "   5) { Id: 48, P: 20, Q: 25, R: 22 }\r\n",
       "   6) { Id: 57, P: 25, Q: 28, R: 26 }\r\n",
       "   7) { Id: 61, P: 27, Q: 30, R: 28 }\r\n",
       "   8) { Id: 68, P: 30, Q: 37, R: 33 }\r\n",
       "   9) { Id: 73, P: 37, Q: 44, R: 40 }\r\n",
       "Seq<{i8,i8,i8,i8}>\r\n",
       "   0) { Id: 0, P: 2, Q: 5, R: 3 }\r\n",
       "   1) { Id: 17, P: 5, Q: 12, R: 8 }\r\n",
       "   2) { Id: 23, P: 6, Q: 11, R: 8 }\r\n",
       "   3) { Id: 36, P: 12, Q: 17, R: 14 }\r\n",
       "   4) { Id: 43, P: 17, Q: 20, R: 18 }\r\n",
       "   5) { Id: 48, P: 20, Q: 25, R: 22 }\r\n",
       "   6) { Id: 57, P: 25, Q: 28, R: 26 }\r\n",
       "   7) { Id: 64, P: 28, Q: 31, R: 29 }\r\n",
       "   8) { Id: 67, P: 30, Q: 33, R: 31 }\r\n",
       "   9) { Id: 74, P: 37, Q: 46, R: 41 }\r\n",
       "Seq<{i8,i8,i8,i8}>\r\n",
       "   0) { Id: 0, P: 2, Q: 5, R: 3 }\r\n",
       "   1) { Id: 17, P: 5, Q: 12, R: 8 }\r\n",
       "   2) { Id: 23, P: 6, Q: 11, R: 8 }\r\n",
       "   3) { Id: 36, P: 12, Q: 17, R: 14 }\r\n",
       "   4) { Id: 43, P: 17, Q: 20, R: 18 }\r\n",
       "   5) { Id: 48, P: 20, Q: 25, R: 22 }\r\n",
       "   6) { Id: 57, P: 25, Q: 28, R: 26 }\r\n",
       "   7) { Id: 61, P: 27, Q: 30, R: 28 }\r\n",
       "   8) { Id: 67, P: 30, Q: 33, R: 31 }\r\n",
       "   9) { Id: 72, P: 33, Q: 46, R: 39 }\r\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Best1.MatchTable;\n",
    "Best2.MatchTable;\n",
    "Best3.MatchTable;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "9e776f53-1756-4d37-8f74-0902eb7c29b3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Seq<bool>\r\n",
       "   0) true\r\n",
       "   1) true\r\n",
       "   2) true\r\n",
       "   3) true\r\n",
       "   4) true\r\n",
       "   5) true\r\n",
       "   6) true\r\n",
       "   7) false\r\n",
       "   8) false\r\n",
       "   9) false\r\n",
       "Seq<bool>\r\n",
       "   0) true\r\n",
       "   1) true\r\n",
       "   2) true\r\n",
       "   3) true\r\n",
       "   4) true\r\n",
       "   5) true\r\n",
       "   6) true\r\n",
       "   7) true\r\n",
       "   8) false\r\n",
       "   9) false\r\n",
       "Seq<bool>\r\n",
       "   0) true\r\n",
       "   1) true\r\n",
       "   2) true\r\n",
       "   3) true\r\n",
       "   4) true\r\n",
       "   5) true\r\n",
       "   6) true\r\n",
       "   7) false\r\n",
       "   8) true\r\n",
       "   9) false\r\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Best1.MatchTable = Best2.MatchTable;\n",
    "Best1.MatchTable = Best3.MatchTable;\n",
    "Best2.MatchTable = Best3.MatchTable;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c5e7bcc-7ba4-4d92-8511-a60d8041401c",
   "metadata": {},
   "source": [
    "### Double\n",
    "\n",
    "Now double the length of the protein and optimize."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "640e9810-3efc-4802-afa3-7e8e8f61043e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "100\r\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Double := Folding=>{ Num: Num * 2, HydroPhobic: HydroPhobic ++ (Num + HydroPhobic) };\n",
    "Double.Num;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "16afd554-619c-46d7-9fa1-4fd51dd98752",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Solver: HiGHS\r\n",
       "20\r\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Wall time: 4530.2879ms"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#!time\n",
    "Best1 := Double->Maximize(NumMatches, \"highs\");\n",
    "Best1.NumMatches;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "27b5a644-bc73-4792-a0a3-2e7f2c645f45",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Solver: Gurobi\r\n",
       "20\r\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Wall time: 49301.131ms"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#!time\n",
    "Best2 := Double->Maximize(NumMatches, \"gurobi\");\n",
    "Best2.NumMatches;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "84c74c8d-3fd5-4f30-aa9c-4eccb7c0eedd",
   "metadata": {},
   "outputs": [],
   "source": [
    "// GLPK struggles with this. I haven't let it run to completion.\n",
    "// #!time\n",
    "// Best3 := Double->Maximize(NumMatches, \"glpk\");\n",
    "// Best3.NumMatches;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1583715e-47ce-485e-bcec-7e8ba246a0d6",
   "metadata": {},
   "source": [
    "## Different Protein Input\n",
    "\n",
    "Here is a second module formulation that takes a single parameter named `Protein` as a text value\n",
    "containg `H` and `-` characters. The `Num` and `HydroPhobic` symbols are now constants."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "e1359a86-e5d6-40d1-add6-60b5ca57f23f",
   "metadata": {},
   "outputs": [],
   "source": [
    "Folding2 := module {\n",
    "    param Protein := \"--H-HHH----HH----H--HH---H-HH-HH-H---H------H-H---\";\n",
    "\n",
    "    // The number of amino acids in the protein.\n",
    "    // The positions of the hydro-phobic amino acids.\n",
    "    const Num := Protein.Len;\n",
    "    const HydroPhobic := Range(Num)->TakeIf(Protein[it:*1] = \"H\");\n",
    "    const NumHydroPhobic := HydroPhobic->Count();\n",
    "\n",
    "    // Construct a table of possible matches, including the matching positions (P and Q)\n",
    "    // and the fold location (R). Then augment with an assigned index (Id) and with\n",
    "    // the sequence of bad fold locations (X).\n",
    "    const PossibleMatches :=\n",
    "        CrossJoin(P: HydroPhobic, Q: HydroPhobic,\n",
    "            Q - P >= 3 and (Q - P) mod 2 = 1, // The join predicate.\n",
    "            { P, Q, R: (P + Q - 1) div 2 })   // The result.\n",
    "        // Add the index (Id) and bad fold locations (X).\n",
    "        +>{ Id: #, X: Range(P, Q)->TakeIf(it != R) };\n",
    "\n",
    "    // The number of possible matches.\n",
    "    const NumPossibleMatches := PossibleMatches->Count();\n",
    "\n",
    "    // Define the free variables:\n",
    "    // * A true/false variable for each possible match.\n",
    "    // * A true/false variable for each possible fold location.\n",
    "    var Matches := Tensor.Fill(false, NumPossibleMatches);\n",
    "    var Folds := Tensor.Fill(false, Num);\n",
    "\n",
    "    // We want to maximize the number of matches.\n",
    "    msr NumMatches := Sum(Matches.Values);\n",
    "\n",
    "    // Construct the constraints.\n",
    "    // For a match, there must be a fold at R.\n",
    "    // For a match, there must NOT be a fold at any of the bad fold locations (X).\n",
    "    con MustFold := PossibleMatches->ForEach(Matches[Id] <= Folds[R]);\n",
    "    con CantFold := PossibleMatches->ForEach(Matches[Id] + Folds[X] <= 1);\n",
    "\n",
    "    // Define the table of matches.\n",
    "    let MatchTable := PossibleMatches->TakeIf(Matches[Id])->{ Id, P, Q, R };\n",
    "};"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "19c4972b-bfb4-48bb-9ac7-25801470e3d3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "50\r\n",
       "18\r\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Folding2.Num;\n",
    "Folding2.NumHydroPhobic;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "8f1af3ee-0e60-4265-902b-9e8314c44b22",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Solver: HiGHS\r\n",
       "10\r\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Wall time: 67.3397ms"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#!time\n",
    "Best := Folding2->Maximize(NumMatches, \"highs\");\n",
    "Best.NumMatches;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8152353d-8d20-48fe-8c48-07b9985b6238",
   "metadata": {},
   "source": [
    "## Double\n",
    "\n",
    "Now double the length of the protein and optimize. Note the simplicity of the new\n",
    "`Protein` value in terms of the old."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "6f29909d-915f-40f7-9ae3-7e8e905ddc76",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "100\r\n",
       "--H-HHH----HH----H--HH---H-HH-HH-H---H------H-H-----H-HHH----HH----H--HH---H-HH-HH-H---H------H-H---\r\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Double := Folding2=>{ Protein: Protein & Protein };\n",
    "Double.Num;\n",
    "Double.Protein;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "0fb8afd1-fe94-4923-b10d-40d0e954bf95",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Solver: HiGHS\r\n",
       "20\r\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Wall time: 4879.3951ms"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#!time\n",
    "Best := Double->Maximize(NumMatches, \"highs\");\n",
    "Best.NumMatches;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "504bd58b-4e15-4095-a6dd-2cd536d5c1be",
   "metadata": {},
   "source": [
    "## Double and Reverse\n",
    "\n",
    "This time, concatenate the original Protein text with it's reverse. Does this allow better matching?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "ebb53fa4-c4db-44bd-8805-efeb05611408",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "100\r\n",
       "--H-HHH----HH----H--HH---H-HH-HH-H---H------H-H------H-H------H---H-HH-HH-H---HH--H----HH----HHH-H--\r\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Double := Folding2=>{ Protein: Protein & Protein[::-1] };\n",
    "Double.Num;\n",
    "Double.Protein;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "66887c27-fa91-4209-81d3-8d16f7c28eea",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Solver: HiGHS\r\n",
       "22\r\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Wall time: 4113.9823ms"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#!time\n",
    "Best := Double->Maximize(NumMatches, \"highs\");\n",
    "Best.NumMatches;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55452a0c-094b-43bd-9841-7862e9dc3fa7",
   "metadata": {},
   "source": [
    "Now in the opposite order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "0791429f-f87e-41f1-b02c-cf9f14247110",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "100\r\n",
       "---H-H------H---H-HH-HH-H---HH--H----HH----HHH-H----H-HHH----HH----H--HH---H-HH-HH-H---H------H-H---\r\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Double := Folding2=>{ Protein: Protein[::-1] & Protein };\n",
    "Double.Num;\n",
    "Double.Protein;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "4182bad3-07bb-4d36-bf6d-f4a6ea8fd6e8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Solver: HiGHS\r\n",
       "22\r\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Wall time: 4119.0905ms"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#!time\n",
    "Best := Double->Maximize(NumMatches, \"highs\");\n",
    "Best.NumMatches;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f75ad828-9aac-493f-8ab2-68be36acd37c",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Rexl",
   "language": "Rexl",
   "name": "rexl"
  },
  "language_info": {
   "file_extension": ".rexl",
   "mimetype": "text/x-rexl",
   "name": "Rexl",
   "pygments_lexer": "rexl",
   "version": "1.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}