From 50dcf7d61ccd3fdf2387d3709c3c0d5f1fb62469 Mon Sep 17 00:00:00 2001 From: Joshua Fogg Date: Tue, 16 Apr 2024 17:43:38 +0100 Subject: [PATCH] Added Julian's n = 2 example --- Example2/A02.txt | 3 + Example2/EBV02.txt | 2 + Example2/README.md | 3 + Example2/S02.txt | 3 + Example2/n2-test.ipynb | 379 +++++++++++++++++++++++++++++++++++++++++ 5 files changed, 390 insertions(+) create mode 100644 Example2/A02.txt create mode 100644 Example2/EBV02.txt create mode 100644 Example2/README.md create mode 100644 Example2/S02.txt create mode 100644 Example2/n2-test.ipynb diff --git a/Example2/A02.txt b/Example2/A02.txt new file mode 100644 index 0000000..3ad0144 --- /dev/null +++ b/Example2/A02.txt @@ -0,0 +1,3 @@ +1 1 1 +1 2 0 +2 2 1 diff --git a/Example2/EBV02.txt b/Example2/EBV02.txt new file mode 100644 index 0000000..1191247 --- /dev/null +++ b/Example2/EBV02.txt @@ -0,0 +1,2 @@ +1 +2 diff --git a/Example2/README.md b/Example2/README.md new file mode 100644 index 0000000..1a3bd3e --- /dev/null +++ b/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/). diff --git a/Example2/S02.txt b/Example2/S02.txt new file mode 100644 index 0000000..82f748c --- /dev/null +++ b/Example2/S02.txt @@ -0,0 +1,3 @@ +1 1 0.11111111111111 +1 2 0 +2 2 4.0 diff --git a/Example2/n2-test.ipynb b/Example2/n2-test.ipynb new file mode 100644 index 0000000..c4fb745 --- /dev/null +++ b/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 +}