Added Julian's n = 2 example

Example2/A02.txt

@@ -0,0 +1,3 @@
+1 1 1
+1 2 0
+2 2 1

Example2/EBV02.txt

@@ -0,0 +1,2 @@
+1
+2

Example2/README.md

@@ -0,0 +1,3 @@
+# Example 2
+
+This is a small example for $n = 2$ created by Julian. It uses the same file formats as [Example 1](../Example/).

Example2/S02.txt

@@ -0,0 +1,3 @@
+1 1 0.11111111111111
+1 2 0
+2 2 4.0

Example2/n2-test.ipynb

@@ -0,0 +1,379 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Julian's Small Example\n",
+    "\n",
+    "Julian came up with a small $n = 2$ example (excluding sex data) for testing the problem. That's solved in this notebook, which requires the following imports."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np                  # defines matrix structures\n",
+    "from qpsolvers import solve_qp      # used for quadratic optimization\n",
+    "from time import perf_counter       # fine grained timing \n",
+    "import gurobipy as gp               # Gurobi optimization interface (1)\n",
+    "from gurobipy import GRB            # Gurobi optimization interface (2)\n",
+    "\n",
+    "# want to round rather than truncate when printing\n",
+    "np.set_printoptions(threshold=np.inf)\n",
+    "\n",
+    "# only show numpy output to five decimal places\n",
+    "np.set_printoptions(formatter={'float_kind':\"{:.5f}\".format})"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Problem Variables\n",
+    "\n",
+    "The variables for Julian's $n = 2$ problem are as follows:\n",
+    "$$\n",
+    "    \\Sigma = \\begin{bmatrix} 1 & 0 \\\\ 0 & 1 \\end{bmatrix},\\quad\\\n",
+    "    \\bar{\\mu} = \\begin{bmatrix} 1 \\\\ 2 \\end{bmatrix},\\quad\\\n",
+    "    \\Omega = \\begin{bmatrix} \\frac{1}{9} & 0 \\\\ 0 & 4 \\end{bmatrix}.\n",
+    "$$\n",
+    "These are stored in `A02.txt`, `EBV02.txt`, and `S02.txt` respectively, where the first and last are matrices in coordinate format. Note that this particular problem does not contain any sex data or weight bound data. The problem is read into Python using the following function."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def load_problem(A_filename, E_filename, S_filename, dimension=False):\n",
+    "    \"\"\"\n",
+    "    Function which loads genetic selection problems in from a\n",
+    "    series of three data files.\n",
+    "    \"\"\"\n",
+    "\n",
+    "    def load_symmetric_matrix(filename, dimension):\n",
+    "        \"\"\"\n",
+    "        Since NumPy doesn't have a stock way to load matrices\n",
+    "        stored in coordinate format format, this adds one.\n",
+    "        \"\"\"\n",
+    "\n",
+    "        matrix = np.zeros([dimension, dimension])\n",
+    "\n",
+    "        with open(filename, 'r') as file:\n",
+    "            for line in file:\n",
+    "                i, j, entry = line.split(\" \")\n",
+    "                # data files indexed from 1, not 0\n",
+    "                matrix[int(i)-1, int(j)-1] = entry\n",
+    "                matrix[int(j)-1, int(i)-1] = entry\n",
+    "\n",
+    "        return matrix\n",
+    "\n",
+    "\n",
+    "    # if dimension wasn't supplied, need to find that\n",
+    "    if not dimension:\n",
+    "        # get dimension from EBV, since it's the smallest file\n",
+    "        with open(E_filename, 'r') as file:\n",
+    "            dimension = sum(1 for _ in file)\n",
+    "\n",
+    "    # EBV isn't in coordinate format so can be loaded directly\n",
+    "    E = np.loadtxt(E_filename)  \n",
+    "    # A and S are stored by coordinates so need special loader\n",
+    "    A = load_symmetric_matrix(A_filename, dimension)\n",
+    "    S = load_symmetric_matrix(S_filename, dimension)\n",
+    "\n",
+    "    return A, E, S, dimension"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The standard problem with $\\lambda = 0.5$ can then be solved using Gurobi as follows:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Set parameter Username\n",
+      "Academic license - for non-commercial use only - expires 2025-02-26\n",
+      "Gurobi Optimizer version 11.0.0 build v11.0.0rc2 (linux64 - \"Ubuntu 22.04.4 LTS\")\n",
+      "\n",
+      "CPU model: Intel(R) Core(TM) i5-8350U CPU @ 1.70GHz, instruction set [SSE2|AVX|AVX2]\n",
+      "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
+      "\n",
+      "Optimize a model with 1 rows, 2 columns and 2 nonzeros\n",
+      "Model fingerprint: 0x37f20dc3\n",
+      "Model has 2 quadratic objective terms\n",
+      "Coefficient statistics:\n",
+      "  Matrix range     [1e+00, 1e+00]\n",
+      "  Objective range  [5e-01, 1e+00]\n",
+      "  QObjective range [1e+00, 1e+00]\n",
+      "  Bounds range     [0e+00, 0e+00]\n",
+      "  RHS range        [1e+00, 1e+00]\n",
+      "Presolve time: 0.01s\n",
+      "Presolved: 1 rows, 2 columns, 2 nonzeros\n",
+      "Presolved model has 2 quadratic objective terms\n",
+      "Ordering time: 0.00s\n",
+      "\n",
+      "Barrier statistics:\n",
+      " AA' NZ     : 0.000e+00\n",
+      " Factor NZ  : 1.000e+00\n",
+      " Factor Ops : 1.000e+00 (less than 1 second per iteration)\n",
+      " Threads    : 1\n",
+      "\n",
+      "                  Objective                Residual\n",
+      "Iter       Primal          Dual         Primal    Dual     Compl     Time\n",
+      "   0   9.99499500e+05 -1.00100025e+06  2.00e+03 0.00e+00  1.00e+06     0s\n",
+      "   1  -4.99749639e-01 -9.99502094e+02  1.00e+00 3.41e-13  9.99e+02     0s\n",
+      "   2  -5.00116624e-01 -2.65997154e+02  1.00e-06 6.99e-15  1.33e+02     0s\n",
+      "   3  -5.00816469e-01 -1.01361177e+00  9.33e-10 0.00e+00  2.56e-01     0s\n",
+      "   4  -5.58837435e-01 -5.96696571e-01  3.43e-12 0.00e+00  1.89e-02     0s\n",
+      "   5  -5.62490657e-01 -5.64212774e-01  0.00e+00 0.00e+00  8.61e-04     0s\n",
+      "   6  -5.62500000e-01 -5.62501707e-01  0.00e+00 0.00e+00  8.53e-07     0s\n",
+      "   7  -5.62500000e-01 -5.62500002e-01  0.00e+00 0.00e+00  8.53e-10     0s\n",
+      "\n",
+      "Barrier solved model in 7 iterations and 0.02 seconds (0.00 work units)\n",
+      "Optimal objective -5.62500000e-01\n",
+      "\n",
+      "Standard: [0.25000 0.75000]\n"
+     ]
+    }
+   ],
+   "source": [
+    "sigma, mubar, omega, n = load_problem(\n",
+    "    \"A02.txt\",\n",
+    "    \"EBV02.txt\",\n",
+    "    \"S02.txt\",\n",
+    "    2)\n",
+    "\n",
+    "lam = 0.5\n",
+    "\n",
+    "# initialise standard genetic selection model\n",
+    "model_std = gp.Model(\"n02standard\")\n",
+    "w_std = model_std.addMVar(shape=n, vtype=GRB.CONTINUOUS, name=\"w\")\n",
+    "\n",
+    "# define the objective functions for standard problem\n",
+    "model_std.setObjective(\n",
+    "    0.5*w_std@(sigma@w_std) - lam*w_std.transpose()@mubar,\n",
+    "GRB.MINIMIZE)\n",
+    "\n",
+    "# add sum-to-one constraint\n",
+    "model_std.addConstr(np.ones([1,2]) @ w_std == 1, name=\"sum-to-one\")\n",
+    "\n",
+    "# solve problem with Gurobi\n",
+    "model_std.optimize()\n",
+    "print(f\"Standard: {w_std.X}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "And the robust problem with $\\kappa = 1.0$ can be solved as follows:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Gurobi Optimizer version 11.0.0 build v11.0.0rc2 (linux64 - \"Ubuntu 22.04.4 LTS\")\n",
+      "\n",
+      "CPU model: Intel(R) Core(TM) i5-8350U CPU @ 1.70GHz, instruction set [SSE2|AVX|AVX2]\n",
+      "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
+      "\n",
+      "Optimize a model with 2 rows, 3 columns and 3 nonzeros\n",
+      "Model fingerprint: 0xb83c03fb\n",
+      "Model has 2 quadratic objective terms\n",
+      "Model has 2 quadratic constraints\n",
+      "Coefficient statistics:\n",
+      "  Matrix range     [1e+00, 1e+00]\n",
+      "  QMatrix range    [1e-01, 4e+00]\n",
+      "  Objective range  [5e-01, 1e+00]\n",
+      "  QObjective range [1e+00, 1e+00]\n",
+      "  Bounds range     [0e+00, 0e+00]\n",
+      "  RHS range        [1e+00, 1e+00]\n",
+      "Presolve removed 1 rows and 0 columns\n",
+      "Presolve time: 0.01s\n",
+      "Presolved: 8 rows, 8 columns, 13 nonzeros\n",
+      "Presolved model has 3 second-order cone constraints\n",
+      "Ordering time: 0.00s\n",
+      "\n",
+      "Barrier statistics:\n",
+      " AA' NZ     : 1.900e+01\n",
+      " Factor NZ  : 3.600e+01\n",
+      " Factor Ops : 2.040e+02 (less than 1 second per iteration)\n",
+      " Threads    : 1\n",
+      "\n",
+      "                  Objective                Residual\n",
+      "Iter       Primal          Dual         Primal    Dual     Compl     Time\n",
+      "   0  -6.30136986e-01 -6.30136986e-01  1.61e+00 8.66e-01  7.30e-01     0s\n",
+      "   1  -1.14616869e+00 -1.60124001e+00  4.19e-01 9.52e-07  1.64e-01     0s\n",
+      "   2  -6.51626451e-01 -8.34836457e-01  4.61e-07 1.05e-12  1.67e-02     0s\n",
+      "   3  -6.72758817e-01 -6.83613315e-01  2.69e-08 3.87e-14  9.87e-04     0s\n",
+      "   4  -6.73595613e-01 -6.73729333e-01  2.55e-14 2.00e-15  1.22e-05     0s\n",
+      "   5  -6.73606055e-01 -6.73616513e-01  1.05e-12 9.30e-13  9.51e-07     0s\n",
+      "   6  -6.73611064e-01 -6.73611270e-01  6.32e-11 1.09e-11  1.88e-08     0s\n",
+      "\n",
+      "Barrier solved model in 6 iterations and 0.03 seconds (0.00 work units)\n",
+      "Optimal objective -6.73611064e-01\n",
+      "\n",
+      "Standard: [0.25000 0.75000]\n",
+      "Robust:   [0.41660 0.58340]\n"
+     ]
+    }
+   ],
+   "source": [
+    "# other than kappa, reuse problem variables from the last cell\n",
+    "kappa = 1.0\n",
+    "\n",
+    "# initialise robust genetic selection model\n",
+    "model_rbs = gp.Model(\"n02robust\")\n",
+    "\n",
+    "# initialise w for both models, z for robust model \n",
+    "w_rbs = model_rbs.addMVar(shape=n, vtype=GRB.CONTINUOUS, name=\"w\")\n",
+    "z_rbs = model_rbs.addVar(name=\"z\")\n",
+    "\n",
+    "model_rbs.setObjective(\n",
+    "    0.5*w_rbs@(sigma@w_rbs) - lam*w_rbs.transpose()@mubar - kappa*z_rbs,\n",
+    "GRB.MINIMIZE)\n",
+    "\n",
+    "# add sum-to-half constraints to both models\n",
+    "model_rbs.addConstr(np.ones([1,2]) @ w_std == 1, name=\"sum-to-one\")\n",
+    "# add quadratic uncertainty constraint to the robust model\n",
+    "model_rbs.addConstr(z_rbs**2 <= np.inner(w_rbs, omega@w_rbs), name=\"uncertainty\")\n",
+    "model_rbs.addConstr(z_rbs >= 0, name=\"z positive\")\n",
+    "\n",
+    "# solve problem with Gurobi, print alongside standard solution for comparison\n",
+    "model_rbs.optimize()\n",
+    "print(f\"Standard: {w_std.X}\")\n",
+    "print(f\"Robust:   {w_rbs.X}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Note that if $\\kappa = 0$, the standard solution not _exactly_ recovered, only accurate to four decimal places."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Gurobi Optimizer version 11.0.0 build v11.0.0rc2 (linux64 - \"Ubuntu 22.04.4 LTS\")\n",
+      "\n",
+      "CPU model: Intel(R) Core(TM) i5-8350U CPU @ 1.70GHz, instruction set [SSE2|AVX|AVX2]\n",
+      "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
+      "\n",
+      "Optimize a model with 2 rows, 3 columns and 3 nonzeros\n",
+      "Model fingerprint: 0x48b37c4b\n",
+      "Model has 2 quadratic objective terms\n",
+      "Model has 2 quadratic constraints\n",
+      "Coefficient statistics:\n",
+      "  Matrix range     [1e+00, 1e+00]\n",
+      "  QMatrix range    [1e-01, 4e+00]\n",
+      "  Objective range  [5e-01, 1e+00]\n",
+      "  QObjective range [1e+00, 1e+00]\n",
+      "  Bounds range     [0e+00, 0e+00]\n",
+      "  RHS range        [1e+00, 1e+00]\n",
+      "Presolve removed 1 rows and 0 columns\n",
+      "Presolve time: 0.01s\n",
+      "Presolved: 8 rows, 8 columns, 13 nonzeros\n",
+      "Presolved model has 3 second-order cone constraints\n",
+      "Ordering time: 0.00s\n",
+      "\n",
+      "Barrier statistics:\n",
+      " AA' NZ     : 1.900e+01\n",
+      " Factor NZ  : 3.600e+01\n",
+      " Factor Ops : 2.040e+02 (less than 1 second per iteration)\n",
+      " Threads    : 1\n",
+      "\n",
+      "                  Objective                Residual\n",
+      "Iter       Primal          Dual         Primal    Dual     Compl     Time\n",
+      "   0  -6.30136986e-01 -6.30136986e-01  1.67e+00 4.15e-01  3.76e-01     0s\n",
+      "   1  -7.63463502e-01 -9.90451466e-01  3.34e-01 4.57e-07  6.75e-02     0s\n",
+      "   2  -5.29769678e-01 -6.77289143e-01  3.67e-07 5.02e-13  1.34e-02     0s\n",
+      "   3  -5.47218963e-01 -5.75552518e-01  4.04e-13 3.33e-16  2.58e-03     0s\n",
+      "   4  -5.61928380e-01 -5.63204951e-01  4.66e-15 1.11e-16  1.16e-04     0s\n",
+      "   5  -5.62460343e-01 -5.62501958e-01  9.44e-14 2.22e-16  3.78e-06     0s\n",
+      "   6  -5.62499955e-01 -5.62500028e-01  2.26e-11 6.10e-14  6.60e-09     0s\n",
+      "\n",
+      "Barrier solved model in 6 iterations and 0.03 seconds (0.00 work units)\n",
+      "Optimal objective -5.62499955e-01\n",
+      "\n",
+      "Standard: [0.2500000031 0.7499999969]\n",
+      "Robust:   [0.2499847190 0.7500152810]\n"
+     ]
+    }
+   ],
+   "source": [
+    "# other than kappa, reuse problem variables from the last cell\n",
+    "kappa = 0.0\n",
+    "\n",
+    "# initialise robust genetic selection model\n",
+    "model_rbs = gp.Model(\"n02robust_k0\")\n",
+    "\n",
+    "# initialise w for both models, z for robust model \n",
+    "w_rbs = model_rbs.addMVar(shape=n, vtype=GRB.CONTINUOUS, name=\"w\")\n",
+    "z_rbs = model_rbs.addVar(name=\"z\")\n",
+    "\n",
+    "model_rbs.setObjective(\n",
+    "    0.5*w_rbs@(sigma@w_rbs) - lam*w_rbs.transpose()@mubar - kappa*z_rbs,\n",
+    "GRB.MINIMIZE)\n",
+    "\n",
+    "# add sum-to-half constraints to both models\n",
+    "model_rbs.addConstr(np.ones([1,2]) @ w_std == 1, name=\"sum-to-one\")\n",
+    "# add quadratic uncertainty constraint to the robust model\n",
+    "model_rbs.addConstr(z_rbs**2 <= np.inner(w_rbs, omega@w_rbs), name=\"uncertainty\")\n",
+    "model_rbs.addConstr(z_rbs >= 0, name=\"z positive\")\n",
+    "\n",
+    "# solve problem with Gurobi, print alongside standard solution for comparison\n",
+    "model_rbs.optimize()\n",
+    "np.set_printoptions(formatter={'float_kind':\"{:.10f}\".format})\n",
+    "print(f\"Standard: {w_std.X}\")\n",
+    "print(f\"Robust:   {w_rbs.X}\")"
+   ]
+  }
+ ],
+ "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.10.12"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}