From fc766b6ad46af0ed369504fec596f0c283eab9d6 Mon Sep 17 00:00:00 2001 From: Joshua Fogg Date: Thu, 18 Jul 2024 18:13:38 +0100 Subject: [PATCH] Switched out in-repo examples for test scripts Example use now lives in the wiki at github.com/Foggalong/RobustOCS/wiki --- examples/04/S04.txt | 2 +- examples/04/example.ipynb | 177 ------ examples/04/test.py | 52 ++ examples/1000/README.md | 11 - examples/1000/solution_rob.txt | 1000 ++++++++++++++++++++++++++++++++ examples/1000/solution_std.txt | 1000 ++++++++++++++++++++++++++++++++ examples/1000/test.py | 50 ++ examples/1000/timing.sh | 38 -- examples/50/example.py | 53 -- examples/50/solution_rob.txt | 50 ++ examples/50/solution_std.txt | 50 ++ examples/50/test.py | 50 ++ examples/README.md | 34 +- 13 files changed, 2255 insertions(+), 312 deletions(-) delete mode 100644 examples/04/example.ipynb create mode 100644 examples/04/test.py delete mode 100644 examples/1000/README.md create mode 100644 examples/1000/solution_rob.txt create mode 100644 examples/1000/solution_std.txt create mode 100644 examples/1000/test.py delete mode 100755 examples/1000/timing.sh delete mode 100644 examples/50/example.py create mode 100644 examples/50/solution_rob.txt create mode 100644 examples/50/solution_std.txt create mode 100644 examples/50/test.py diff --git a/examples/04/S04.txt b/examples/04/S04.txt index ae2331b..017b797 100644 --- a/examples/04/S04.txt +++ b/examples/04/S04.txt @@ -1,4 +1,4 @@ 1 1 0.11111111111111 2 2 0.11111111111111 3 3 4.0 -4 4 4.0 \ No newline at end of file +4 4 4.0 diff --git a/examples/04/example.ipynb b/examples/04/example.ipynb deleted file mode 100644 index 59991e3..0000000 --- a/examples/04/example.ipynb +++ /dev/null @@ -1,177 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Julian's Small Example\n", - "\n", - "[Julian](https://www.maths.ed.ac.uk/hall) came up with a small $n = 2$ example (excluding sex data) for testing the theory behind solving the problem. This example can easily be extended and then solved using the RobustOCS module. This is a quick demonstration of using the module within an interactive notebook environment." - ] - }, - { - "cell_type": "code", - "execution_count": 1, - "metadata": {}, - "outputs": [], - "source": [ - "import robustocs" - ] - }, - { - "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", - "This can easily be extended to include sex data by reflecting those variable across both sires and dams, filling in the remainder with zeros. \n", - "$$\n", - " \\Sigma = \\begin{bmatrix}\n", - " 1 & 0 & 0 & 0 \\\\\n", - " 0 & 1 & 0 & 0 \\\\\n", - " 0 & 0 & 1 & 0 \\\\\n", - " 0 & 0 & 0 & 1\n", - " \\end{bmatrix},\\quad\\\n", - " \\bar{\\mu} = \\begin{bmatrix} 1 \\\\ 1 \\\\ 2 \\\\ 2 \\end{bmatrix},\\quad\\\n", - " \\Omega = \\begin{bmatrix}\n", - " \\frac{1}{9} & 0 & 0 & 0 \\\\\n", - " 0 & \\frac{1}{9} & 0 & 0 \\\\\n", - " 0 & 0 & 4 & 0 \\\\\n", - " 0 & 0 & 0 & 4\n", - " \\end{bmatrix},\\quad\\\n", - " \\mathcal{S} = \\lbrace 1, 3 \\rbrace,\\quad\\\n", - " \\mathcal{D} = \\lbrace 2, 4 \\rbrace.\n", - "$$\n", - "The matrix variables are stored in [`A04.txt`](A04.txt), [`EBV04.txt`](EBV04.txt), and [`S04.txt`](S04.txt) respectively, where the first and last are matrices in coordinate format. These are then loaded into Python using `load_problem`." - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": {}, - "outputs": [], - "source": [ - "# key problem variables loaded from standard format txt files\n", - "sigma, mubar, omega, n = robustocs.load_problem(\"A04.txt\", \"EBV04.txt\", \"S04.txt\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Since we have constructed the problem to have alternating sex data, we may approach it as we did in [`example.py`](../50/example.py) using range iterates." - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": {}, - "outputs": [], - "source": [ - "sires = range(0, n, 2)\n", - "dams = range(1, n, 2)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Gurobi\n", - "\n", - "AlphaRGS defines functions for solving the standard and robust genetic selection problems using the [gurobipy](https://pypi.org/project/gurobipy/) Python interface to Gurobi. Below these are used to solve the above problem for $\\lambda = 0.5, \\kappa = 1$: " - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Set parameter Username\n", - "Academic license - for non-commercial use only - expires 2025-02-26\n", - "i w_std w_rbs\n", - "1 0.00000 0.38200\n", - "2 0.00000 0.38200\n", - "3 0.50000 0.11800\n", - "4 0.50000 0.11800\n", - "\n", - "w_std objective: 1.87500\n", - "w_rbs objective: 0.77684 (z = 0.37924)\n", - "Maximum change: 0.38200\n", - "Average change: 0.38200\n", - "Minimum change: 0.38200\n" - ] - } - ], - "source": [ - "lam = 0.5\n", - "kap = 1\n", - "\n", - "# computes the standard and robust genetic selection solutions\n", - "w_std, obj_std = robustocs.gurobi_standard_genetics(sigma, mubar, sires, dams, lam, n)\n", - "w_rbs, z_rbs, obj_rbs = robustocs.gurobi_robust_genetics(sigma, mubar, omega, sires, dams, lam, kap, n)\n", - "\n", - "# print a comparison of the two solutions\n", - "robustocs.print_compare_solutions(w_std, w_rbs, obj_std, obj_rbs, z2=z_rbs, name1=\"w_std\", name2=\"w_rbs\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Finally, we can use `check_uncertainty_constraint` to explore how close our $z\\geq \\sqrt{w^{T}\\Omega w}$ constraint came to equality." - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "\n", - " z: 0.37923871642022844\n", - "w'*Ω*w: 0.3792386953366983\n", - " Diff: 2.1083530143961582e-08\n" - ] - } - ], - "source": [ - "if not robustocs.check_uncertainty_constraint(z_rbs, w_rbs, omega, debug=True):\n", - " raise ValueError" - ] - } - ], - "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 -} diff --git a/examples/04/test.py b/examples/04/test.py new file mode 100644 index 0000000..a48d5db --- /dev/null +++ b/examples/04/test.py @@ -0,0 +1,52 @@ +#!/usr/bin/env python3 + +import numpy as np +import robustocs as rocs + +# SETUP +# ----- + +# load in the problem variables +sigma, mubar, omega, n = rocs.load_problem( + "A04.txt", + "EBV04.txt", + "S04.txt", + issparse=True +) +sires = range(0, n, 2) +dams = range(1, n, 2) +lam = 0.5 +kap = 1 + +# true solution to the standard and robust problems +true_std = np.array([0, 0, 0.5, 0.5]) +true_rob = np.array([0.382, 0.382, 0.118, 0.118]) + +# tolerance for comparisons; NOTE this test uses a much stricter tolerance +# for the non-robust problem since both solvers should compute it exactly. +tol_std = 1e-7 +tol_rob = 1e-3 + + +# TESTS +# ----- + +w, obj = rocs.highs_standard_genetics(sigma, mubar, sires, dams, lam, n) +assert ((w - true_std) < tol_std).all(), "QP in HiGHS was incorrect" + +w, obj = rocs.gurobi_standard_genetics(sigma, mubar, sires, dams, lam, n) +assert ((w - true_std) < tol_std).all(), "QP in Gurobi was incorrect" + +w, z, obj = rocs.gurobi_robust_genetics( + sigma, mubar, omega, sires, dams, lam, kap, n) +assert ((w - true_rob) < tol_rob).all(), "conic in Gurobi was incorrect" + +w, z, obj = rocs.gurobi_robust_genetics_sqp( + sigma, mubar, omega, sires, dams, lam, kap, n) +assert ((w - true_rob) < tol_rob).all(), "SQP in Gurobi was incorrect" + +w, z, obj = rocs.highs_robust_genetics_sqp( + sigma, mubar, omega, sires, dams, lam, kap, n) +assert ((w - true_rob) < tol_rob).all(), "conic in HiGHS was incorrect" + +print("Success!") diff --git a/examples/1000/README.md b/examples/1000/README.md deleted file mode 100644 index c45b3d3..0000000 --- a/examples/1000/README.md +++ /dev/null @@ -1,11 +0,0 @@ -# Timings - -This is how long to solve the $n = 1000$ example with each method at various commits. Timings in seconds (all shown to three significant figures) are the average of three runs, repeated five times with the best average taken. The first two solvers are standard non-robust optimization, and the last three are all robust optimization (conic or SQP). - -| Commit | Gurobi (standard) | HiGHS (standard) | Gurobi (conic) | Gurobi (SQP) | HiGHS (SQP) | Total | -| :-----: | ----------------: | ---------------: | -------------: | -----------: | ----------: | ----: | -| 7d866d8 | 0.688 | 0.197 | 2.690 | 24.400 | 1.700 | 492.8 | -| 90a6040 | 0.690 | 0.204 | 2.770 | 24.600 | 1.700 | 499.9 | -| 92ae275 | 0.676 | 0.204 | 2.750 | 24.600 | 1.680 | 507.3 | -| 48e645f | 0.675 | 0.205 | 2.750 | 24.600 | 1.680 | 508.7 | -| 3433097 | 0.676 | 0.204 | 2.760 | 24.400 | 1.700 | 508.1 | diff --git a/examples/1000/solution_rob.txt b/examples/1000/solution_rob.txt new file mode 100644 index 0000000..7aec2bd --- /dev/null +++ b/examples/1000/solution_rob.txt @@ -0,0 +1,1000 @@ +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0179 +0 +0.0051 +0 +0 +0 +0 +0 +0 +0 +0.0157 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0053 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0313 +0 +0 +0 +0 +0.0362 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0517 +0 +0 +0 +0 +0.0545 +0 +0 +0 +0 +0 +0.0116 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0148 +0 +0 +0 +0.0034 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0171 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0125 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0213 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0004 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0093 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0082 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0034 +0.0197 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0333 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0030 +0 +0 +0 +0 +0 +0 +0 +0.0343 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0149 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0261 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0101 +0 +0.0334 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0009 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0211 +0 +0 +0 +0 +0.0446 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0778 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0043 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0073 +0.0681 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0128 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0224 +0 +0.0768 +0 +0.0002 +0 +0.1501 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0161 +0 +0 +0.0032 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 diff --git a/examples/1000/solution_std.txt b/examples/1000/solution_std.txt new file mode 100644 index 0000000..b15ac84 --- /dev/null +++ b/examples/1000/solution_std.txt @@ -0,0 +1,1000 @@ +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0253 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.1357 +0 +0 +0 +0 +0.0424 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0706 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0810 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0195 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.1679 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0522 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0890 +0 +0 +0 +0.3164 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 diff --git a/examples/1000/test.py b/examples/1000/test.py new file mode 100644 index 0000000..fd26317 --- /dev/null +++ b/examples/1000/test.py @@ -0,0 +1,50 @@ +#!/usr/bin/env python3 + +import numpy as np +import robustocs as rocs + +# SETUP +# ----- + +# load in the problem variables +sigma, mubar, omega, n = rocs.load_problem( + "A1000.txt", + "EBV1000.txt", + "S1000.txt", + issparse=True +) +sires = range(0, n, 2) +dams = range(1, n, 2) +lam = 0.5 +kap = 1 + +# load in true solution to the standard and robust problems +true_std = np.loadtxt('solution_std.txt') +true_rob = np.loadtxt('solution_rob.txt') + +# tolerance for comparisons +tol = 1e-3 + + +# TESTS +# ----- + +w, obj = rocs.highs_standard_genetics(sigma, mubar, sires, dams, lam, n) +assert ((w - true_std) < tol).all(), "QP in HiGHS was incorrect" + +w, obj = rocs.gurobi_standard_genetics(sigma, mubar, sires, dams, lam, n) +assert ((w - true_std) < tol).all(), "QP in Gurobi was incorrect" + +w, z, obj = rocs.gurobi_robust_genetics( + sigma, mubar, omega, sires, dams, lam, kap, n) +assert ((w - true_rob) < tol).all(), "conic in Gurobi was incorrect" + +w, z, obj = rocs.gurobi_robust_genetics_sqp( + sigma, mubar, omega, sires, dams, lam, kap, n) +assert ((w - true_rob) < tol).all(), "SQP in Gurobi was incorrect" + +w, z, obj = rocs.highs_robust_genetics_sqp( + sigma, mubar, omega, sires, dams, lam, kap, n) +assert ((w - true_rob) < tol).all(), "conic in HiGHS was incorrect" + +print("Success!") diff --git a/examples/1000/timing.sh b/examples/1000/timing.sh deleted file mode 100755 index f5b6ae0..0000000 --- a/examples/1000/timing.sh +++ /dev/null @@ -1,38 +0,0 @@ -#!/usr/bin/env bash - -LOOPS=3 # run the command this many times and take the average -REPS=5 # take the best average time from this many averages - -# code block which sets up the problem, same across all three -SETUP=' -import sys; sys.path.insert(1, "../../") -import robustocs as rocs - -sigma, mubar, omega, n = rocs.load_problem( - "A1000.txt", - "EBV1000.txt", - "S1000.txt", - issparse=True -) - -sires = range(0, n, 2) -dams = range(1, n, 2) - -lam = 0.5 -kap = 1' - -# array of solver commands to time -SOLVERS=( - 'rocs.gurobi_standard_genetics(sigma, mubar, sires, dams, lam, n)' - 'rocs.highs_standard_genetics(sigma, mubar, sires, dams, lam, n)' - 'rocs.gurobi_robust_genetics(sigma, mubar, omega, sires, dams, lam, kap, n)' - 'rocs.gurobi_robust_genetics_sqp(sigma, mubar, omega, sires, dams, lam, kap, n)' - 'rocs.highs_robust_genetics_sqp(sigma, mubar, omega, sires, dams, lam, kap, n)' -) - -# run timing code for each solver command -for f in "${SOLVERS[@]}"; do - echo -e "\033[1mTesting ${f}\033[0m" # print command in bold - python3 -m timeit -v -s "$SETUP" -n $LOOPS -r $REPS -u "sec" "$f" - echo -e "\n" # padding for readability -done diff --git a/examples/50/example.py b/examples/50/example.py deleted file mode 100644 index b93a836..0000000 --- a/examples/50/example.py +++ /dev/null @@ -1,53 +0,0 @@ -#!/usr/bin/env python3 - -import sys; sys.path.insert(1, '../../') # HACK pre-module workaround -import robustocs -import time - -# key problem variables loaded from standard format txt files -sigma, mubar, omega, n = robustocs.load_problem( - "A50.txt", - "EBV50.txt", - "S50.txt", - issparse=True -) - -# NOTE this trick of handling sex data is specific to the initial simulation -# data which is structured so that candidates alternate between sires (which -# have even indices) and dams (which have odd indices). -sires = range(0, n, 2) -dams = range(1, n, 2) - -lam = 0.5 -kap = 1 - - -# Gurobi using direct optimization -t0 = time.time() -w_grb, z_grb, obj_grb = robustocs.gurobi_robust_genetics( - sigma, mubar, omega, sires, dams, lam, kap, n) -t1 = time.time() -print(f"Gurobi took {t1-t0:.5f} seconds (direct)") - -# Gurobi using sequential quadratic programming -t0 = time.time() -w_grb, z_grb, obj_grb = robustocs.gurobi_robust_genetics_sqp( - sigma, mubar, omega, sires, dams, lam, kap, n) -t1 = time.time() -print(f"Gurobi took {t1-t0:.5f} seconds (SQP)") - -# HiGHS using sequential quadratic programming -t0 = time.time() -w_hig, z_hig, obj_hig = robustocs.highs_robust_genetics_sqp( - sigma, mubar, omega, sires, dams, lam, kap, n) -t1 = time.time() -print(f"HiGHS took {t1-t0:.5f} seconds (SQP)") - - -print("\nSQP Methods:") -robustocs.print_compare_solutions( - w_grb, w_hig, obj_grb, obj_hig, z1=z_grb, z2=z_hig, - name1="Gurobi", name2="HiGHS ", tol=1e-7 -) - -print("\nDone!") diff --git a/examples/50/solution_rob.txt b/examples/50/solution_rob.txt new file mode 100644 index 0000000..62cba9c --- /dev/null +++ b/examples/50/solution_rob.txt @@ -0,0 +1,50 @@ +0 +0.0781 +0.1247 +0 +0 +0 +0 +0.0808 +0 +0.0353 +0 +0 +0 +0 +0.1059 +0.0671 +0 +0 +0 +0 +0 +0 +0 +0 +0.0340 +0 +0 +0 +0 +0 +0 +0.0446 +0.0934 +0 +0 +0.0104 +0.1420 +0.0791 +0 +0 +0 +0.0814 +0 +0 +0 +0 +0 +0.0232 +0 +0 diff --git a/examples/50/solution_std.txt b/examples/50/solution_std.txt new file mode 100644 index 0000000..b34ed95 --- /dev/null +++ b/examples/50/solution_std.txt @@ -0,0 +1,50 @@ +0 +0.1194 +0.1800 +0 +0 +0 +0 +0.1517 +0 +0 +0 +0 +0 +0 +0.0816 +0.0555 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0.0233 +0 +0 +0 +0.2151 +0.0773 +0 +0 +0 +0.0961 +0 +0 +0 +0 +0 +0 +0 +0 diff --git a/examples/50/test.py b/examples/50/test.py new file mode 100644 index 0000000..c60c700 --- /dev/null +++ b/examples/50/test.py @@ -0,0 +1,50 @@ +#!/usr/bin/env python3 + +import numpy as np +import robustocs as rocs + +# SETUP +# ----- + +# load in the problem variables +sigma, mubar, omega, n = rocs.load_problem( + "A50.txt", + "EBV50.txt", + "S50.txt", + issparse=True +) +sires = range(0, n, 2) +dams = range(1, n, 2) +lam = 0.5 +kap = 1 + +# load in true solution to the standard and robust problems +true_std = np.loadtxt('solution_std.txt') +true_rob = np.loadtxt('solution_rob.txt') + +# tolerance for comparisons +tol = 1e-3 + + +# TESTS +# ----- + +w, obj = rocs.highs_standard_genetics(sigma, mubar, sires, dams, lam, n) +assert ((w - true_std) < tol).all(), "QP in HiGHS was incorrect" + +w, obj = rocs.gurobi_standard_genetics(sigma, mubar, sires, dams, lam, n) +assert ((w - true_std) < tol).all(), "QP in Gurobi was incorrect" + +w, z, obj = rocs.gurobi_robust_genetics( + sigma, mubar, omega, sires, dams, lam, kap, n) +assert ((w - true_rob) < tol).all(), "conic in Gurobi was incorrect" + +w, z, obj = rocs.gurobi_robust_genetics_sqp( + sigma, mubar, omega, sires, dams, lam, kap, n) +assert ((w - true_rob) < tol).all(), "SQP in Gurobi was incorrect" + +w, z, obj = rocs.highs_robust_genetics_sqp( + sigma, mubar, omega, sires, dams, lam, kap, n) +assert ((w - true_rob) < tol).all(), "conic in HiGHS was incorrect" + +print("Success!") diff --git a/examples/README.md b/examples/README.md index c64685c..8dbb288 100644 --- a/examples/README.md +++ b/examples/README.md @@ -1,33 +1,3 @@ -# Examples +# Example Data -- [04/](04/) A scaled up version of the small $n = 2$ example created by Julian, replicating the same problem across sires and dams. -- [50/](50/) The 50 youngest individuals from a generated cohort of 12,000. Implications of this non-random example are unknown. - -## File Description - -In the original simulation, we had 12k individuals and 1k samples of their estimated breeding values (EBV) – our selection criterion. - -Based on that simulated data: - -- I constructed the A-matrix (12k x 12k), this is measure of co-ancestry between individuals based on the pedigree data. This is the standard used in optimal contributions, but it can be constructed in different ways (we can later test the DNA method vs pedigree method, etc.). [A50] -- I constructed the S-matrix (12k x 12k), which is the covariance matrix between 1k EBV samples for 12k individuals. [S50] -- The EBV vector (12k), which is posterior mean over 1k samples of EBV. [EBV50] - -The matrices are saved in row, column, value format like we talked about before. - -_Note: this is an unedited version of Gregor's description of the files._ - -## Times - -To give an idea of how the methods available scale with dimension, the table below times each for four example problems of increasing size. - -| $n$ | Gurobi (standard) | HiGHS (standard) | Gurobi (conic) | Gurobi (SQP) | HiGHS (SQP) | -| -----: | ----------------: | ---------------: | -------------: | -----------: | ----------: | -| 4 | 2.95e-3 | 4.78e-4 | 4.71e-3 | 1.74e-2 | 5.98e-3 | -| 50 | 4.46e-3 | 1.02e-3 | 1.02e-2 | 5.52e-2 | 1.84e-2 | -| 1000 | 6.76e-1 | 2.04e-1 | 2.75e+0 | 2.64e+1 | 1.68e+0 | -| 10000¹ | 8.63e+1 | 2.58e+1 | DNF² | 1.56e+3 | 1.06e+2 | - -_1: This repository doesn't contain the data files due to storage limitations. Contact the authors for access._ - -_2: Gurobi crashed without displaying an error message when attempting to solve using conic programming for the largest problem._ +These are datafiles for example problems, used in the [repository demo](../example.py), in the [GitHub wiki](https://github.com/Foggalong/RobustOCS/wiki), and for testing. A description of the file formats can be [found here](https://github.com/Foggalong/RobustOCS/wiki/File-Formats).