2020_01_12a_nn_least_sq_optimization_with_pyomo_and_ipopt.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To install pyomo and ipopt with conda, execute this in your command line\n",
    "\n",
    "    conda install -c conda-forge pyomo ipopt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal value of x: (0.5, 0.5)\n"
     ]
    }
   ],
   "source": [
    "# here is a code block that will confirm pyomo and ipopt are installed and working\n",
    "\n",
    "\n",
    "import pyomo.environ\n",
    "from pyomo.core import *\n",
    "from pyomo.opt import SolverFactory\n",
    "\n",
    "m = ConcreteModel()\n",
    "m.x = Var([1,2], within=NonNegativeReals)\n",
    "m.objective = Objective(expr=m.x[1]**2 + m.x[2]**2)\n",
    "m.constraint = Constraint(expr=m.x[1] + m.x[2] == 1)\n",
    "    \n",
    "solver = SolverFactory('ipopt')\n",
    "results = solver.solve(m)\n",
    "\n",
    "print('Optimal value of x:', (value(m.x[1]), value(m.x[2])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Too cool!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[13.333333336604015, 13.333333336604015, 13.333333336604015, 0.0]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def nonnegative_optimize(imprecise_counts, control_total):\n",
    "    \"\"\"optimize the imprecise counts so that they sum to\n",
    "    the control total and are non-negative\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    imprecise_counts : list-like of floats\n",
    "    control_total : float\n",
    "    \n",
    "    Results\n",
    "    -------\n",
    "    returns optimized_counts, which are close to imprecise counts,\n",
    "    but not negative, and match control total in aggregate\n",
    "    \"\"\"\n",
    "    imprecise_counts = list(imprecise_counts)\n",
    "    \n",
    "    model = ConcreteModel()\n",
    "    model.I = range(len(imprecise_counts))\n",
    "    model.x = Var(model.I, within=NonNegativeReals)\n",
    "    model.objective = Objective(\n",
    "        expr=sum((model.x[i] - imprecise_counts[i])**2 for i in model.I))\n",
    "    model.constraint = Constraint(\n",
    "        expr=summation(model.x) == control_total)\n",
    "    \n",
    "    solver = SolverFactory('ipopt')\n",
    "    results = solver.solve(model, options={'acceptable_tol':1e-4}, tee=False)\n",
    "    optimized_counts = [value(model.x[i]) for i in model.I]\n",
    "        \n",
    "    return optimized_counts\n",
    "nonnegative_optimize([10, 10, 10, -10], 40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np, matplotlib.pyplot as plt, pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 48 ms, sys: 4 ms, total: 52 ms\n",
      "Wall time: 128 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(array([34.,  6., 24., 12., 12.,  6.,  4.,  1.,  0.,  1.]),\n",
       " array([0.        , 0.4101644 , 0.82032881, 1.23049321, 1.64065761,\n",
       "        2.05082202, 2.46098642, 2.87115082, 3.28131523, 3.69147963,\n",
       "        4.10164403]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAM90lEQVR4nO3df6jd9X3H8edrSVplFlQ8lWBktxTplIJx3GVCYOtSO1Id08IKEyYOhHSgoCDb0v7TFjZIYdX9MwrpdAbm7EQtinY/go2IUOxubJrG3RY7ybbUYK44p/nHkfjeH/ebkd6cm3PuPefccz/N8wGHe873fM/9vvPVPPnyvd9vbqoKSVJ7fmnaA0iSVseAS1KjDLgkNcqAS1KjDLgkNWrjWm7siiuuqJmZmbXcpCQ17+DBg29VVW/p8oEBT3IR8CLw4W79J6rqy0keAX4L+J9u1T+qqkPn+14zMzPMzc2tdHZJuqAl+Y9+y4c5An8f2FFVJ5NsAl5K8o/de39SVU+Ma0hJ0vAGBrwW7/Q52b3c1D28+0eSpmyoH2Im2ZDkEHAC2F9VL3dv/UWSw0keTPLhiU0pSTrHUAGvqtNVtRXYAmxL8kngi8CvAr8OXA78Wb/PJtmVZC7J3MLCwpjGliSt6DLCqnoHeAHYWVXHa9H7wN8C25b5zN6qmq2q2V7vnB+iSpJWaWDAk/SSXNo9vxi4Cfhxks3dsgC3AUcmOagk6ecNcxXKZmBfkg0sBv/xqno2yXeT9IAAh4A/nuCckqQlhrkK5TBwQ5/lOyYykSRpKN5KL0mNWtNb6Ucxs/u5qW376J5bprZtSVqOR+CS1CgDLkmNMuCS1CgDLkmNMuCS1CgDLkmNMuCS1CgDLkmNMuCS1CgDLkmNMuCS1CgDLkmNMuCS1CgDLkmNMuCS1CgDLkmNMuCS1CgDLkmNMuCS1CgDLkmNGhjwJBcl+X6SHyZ5NclXu+UfS/JykteS/EOSD01+XEnSGcMcgb8P7Kiq64GtwM4kNwJfAx6sqmuA/wbumtyYkqSlBga8Fp3sXm7qHgXsAJ7olu8DbpvIhJKkvoY6B55kQ5JDwAlgP/DvwDtVdapb5Rhw1TKf3ZVkLsncwsLCOGaWJDFkwKvqdFVtBbYA24Br+622zGf3VtVsVc32er3VTypJ+jkrugqlqt4BXgBuBC5NsrF7awvwxnhHkySdzzBXofSSXNo9vxi4CZgHDgC/3612J/D0pIaUJJ1r4+BV2AzsS7KBxeA/XlXPJvk34FtJ/hz4AfDQBOeUJC0xMOBVdRi4oc/y11k8Hy5JmgLvxJSkRhlwSWqUAZekRhlwSWqUAZekRhlwSWqUAZekRhlwSWqUAZekRhlwSWqUAZekRhlwSWqUAZekRhlwSWqUAZekRhlwSWqUAZekRhlwSWqUAZekRhlwSWqUAZekRg0MeJKrkxxIMp/k1ST3dsu/kuRnSQ51j5snP64k6YyNQ6xzCri/ql5J8hHgYJL93XsPVtVfTm48SdJyBga8qo4Dx7vn7yWZB66a9GCSpPNb0TnwJDPADcDL3aJ7khxO8nCSy8Y8myTpPIYOeJJLgCeB+6rqXeAbwMeBrSweoX99mc/tSjKXZG5hYWEMI0uSYMiAJ9nEYrwfraqnAKrqzao6XVUfAN8EtvX7bFXtrarZqprt9XrjmluSLnjDXIUS4CFgvqoeOGv55rNW+xxwZPzjSZKWM8xVKNuBO4AfJTnULfsScHuSrUABR4EvTGRCSVJfw1yF8hKQPm99Z/zjSJKG5Z2YktSoYU6haEpmdj83le0e3XPLVLYraWU8ApekRhlwSWqUAZekRhlwSWqUAZekRhlwSWqUAZekRhlwSWqUAZekRhlwSWqUAZekRhlwSWqUAZekRhlwSWqUAZekRhlwSWqUAZekRhlwSWqUAZekRhlwSWrUwIAnuTrJgSTzSV5Ncm+3/PIk+5O81n29bPLjSpLOGOYI/BRwf1VdC9wI3J3kOmA38HxVXQM8372WJK2RgQGvquNV9Ur3/D1gHrgKuBXY1622D7htUkNKks61onPgSWaAG4CXgSur6jgsRh746DKf2ZVkLsncwsLCaNNKkv7f0AFPcgnwJHBfVb077Oeqam9VzVbVbK/XW82MkqQ+hgp4kk0sxvvRqnqqW/xmks3d+5uBE5MZUZLUzzBXoQR4CJivqgfOeusZ4M7u+Z3A0+MfT5K0nI1DrLMduAP4UZJD3bIvAXuAx5PcBfwn8PnJjChJ6mdgwKvqJSDLvP3p8Y4jSRqWd2JKUqMMuCQ1yoBLUqMMuCQ1yoBLUqMMuCQ1yoBLUqMMuCQ1yoBLUqMMuCQ1yoBLUqMMuCQ1yoBLUqMMuCQ1yoBLUqMMuCQ1yoBLUqMMuCQ1yoBLUqMMuCQ1yoBLUqMGBjzJw0lOJDly1rKvJPlZkkPd4+bJjilJWmqYI/BHgJ19lj9YVVu7x3fGO5YkaZCBAa+qF4G312AWSdIKjHIO/J4kh7tTLJctt1KSXUnmkswtLCyMsDlJ0tlWG/BvAB8HtgLHga8vt2JV7a2q2aqa7fV6q9ycJGmpVQW8qt6sqtNV9QHwTWDbeMeSJA2yqoAn2XzWy88BR5ZbV5I0GRsHrZDkMeBTwBVJjgFfBj6VZCtQwFHgCxOcUZLUx8CAV9XtfRY/NIFZJEkr4J2YktQoAy5JjTLgktQoAy5JjTLgktQoAy5JjTLgktQoAy5JjTLgktQoAy5JjTLgktQoAy5JjTLgktSogf8aoS48M7ufm/YIa+7onlumPYK0Yh6BS1KjDLgkNcqAS1KjDLgkNcqAS1KjDLgkNcqAS1KjDLgkNWpgwJM8nOREkiNnLbs8yf4kr3VfL5vsmJKkpYY5An8E2Llk2W7g+aq6Bni+ey1JWkMDA15VLwJvL1l8K7Cve74PuG3Mc0mSBljtOfArq+o4QPf1o8utmGRXkrkkcwsLC6vcnCRpqYn/ELOq9lbVbFXN9nq9SW9Oki4Yqw34m0k2A3RfT4xvJEnSMFYb8GeAO7vndwJPj2ccSdKwhrmM8DHge8AnkhxLchewB/hMkteAz3SvJUlraOAvdKiq25d569NjnkWStALeiSlJjTLgktQoAy5JjTLgktQoAy5JjTLgktQoAy5JjTLgktQoAy5JjTLgktQoAy5JjTLgktQoAy5JjTLgktQoAy5JjTLgktQoAy5JjTLgktSogb9STTCz+7lpj6AJm+Z/46N7bpnattU2j8AlqVEGXJIaNdIplCRHgfeA08Cpqpodx1CSpMHGcQ78t6vqrTF8H0nSCngKRZIaNWrAC/iXJAeT7Oq3QpJdSeaSzC0sLIy4OUnSGaMGfHtV/RrwWeDuJL+5dIWq2ltVs1U12+v1RtycJOmMkQJeVW90X08A3wa2jWMoSdJgqw54kl9O8pEzz4HfAY6MazBJ0vmNchXKlcC3k5z5Pn9fVf80lqkkSQOtOuBV9Tpw/RhnkSStgP8WijRl0/p3WPw3WNrndeCS1CgDLkmNMuCS1CgDLkmNMuCS1CgDLkmNMuCS1CgDLkmNMuCS1CgDLkmNMuCS1CgDLkmNMuCS1CgDLkmNMuCS1CgDLkmNMuCS1CgDLkmN8leqSReoaf0qN5jer3P7RfszewQuSY0y4JLUqJECnmRnkp8k+WmS3eMaSpI02KoDnmQD8NfAZ4HrgNuTXDeuwSRJ5zfKEfg24KdV9XpV/S/wLeDW8YwlSRpklKtQrgL+66zXx4DfWLpSkl3Aru7lySQ/WeX2rgDeWuVnLyTup+G5r4Yz9v2Ur43zu60b591PI/6Zf6XfwlECnj7L6pwFVXuBvSNsZ3FjyVxVzY76fX7RuZ+G574ajvtpONPYT6OcQjkGXH3W6y3AG6ONI0ka1igB/1fgmiQfS/Ih4A+AZ8YzliRpkFWfQqmqU0nuAf4Z2AA8XFWvjm2yc418GuYC4X4anvtqOO6n4az5fkrVOaetJUkN8E5MSWqUAZekRjURcG/ZHyzJw0lOJDky7VnWsyRXJzmQZD7Jq0nunfZM61GSi5J8P8kPu/301WnPtJ4l2ZDkB0meXcvtrvuAe8v+0B4Bdk57iAacAu6vqmuBG4G7/f+pr/eBHVV1PbAV2JnkxinPtJ7dC8yv9UbXfcDxlv2hVNWLwNvTnmO9q6rjVfVK9/w9Fv/SXTXdqdafWnSye7mpe3jFQx9JtgC3AH+z1ttuIeD9btn3L5xGlmQGuAF4ebqTrE/daYFDwAlgf1W5n/r7K+BPgQ/WesMtBHyoW/allUhyCfAkcF9VvTvtedajqjpdVVtZvMt6W5JPTnum9SbJ7wInqurgNLbfQsC9ZV9jlWQTi/F+tKqemvY8611VvQO8gD9j6Wc78HtJjrJ4endHkr9bq423EHBv2dfYJAnwEDBfVQ9Me571KkkvyaXd84uBm4AfT3eq9aeqvlhVW6pqhsU2fbeq/nCttr/uA15Vp4Azt+zPA49P+Jb9JiV5DPge8Ikkx5LcNe2Z1qntwB0sHikd6h43T3uodWgzcCDJYRYPovZX1ZpeIqfBvJVekhq17o/AJUn9GXBJapQBl6RGGXBJapQBl6RGGXBJapQBl6RG/R8VQRSczh23WgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.random.seed(12345)\n",
    "I = 100\n",
    "%time plt.hist(nonnegative_optimize(np.random.normal(1,1,size=I), I))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 108 ms, sys: 20 ms, total: 128 ms\n",
      "Wall time: 402 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(array([347., 185., 184., 135.,  83.,  40.,  21.,   3.,   1.,   1.]),\n",
       " array([0.        , 0.48648472, 0.97296944, 1.45945415, 1.94593887,\n",
       "        2.43242359, 2.91890831, 3.40539302, 3.89187774, 4.37836246,\n",
       "        4.86484718]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAQIklEQVR4nO3df6xfdX3H8edrhaFBN2BcSNc2u8R1mz8SC7lDEpLFgVMEs2IyFkimxJDUJZhg5n4U/1GTkWCisphsJHUwy+bERjQ0wpwdYgzJAG+xVEpldtrJtQ29DkGIGQv1vT/u6byUb/v93vu93/vlfu7zkZx8z/mczznnfdL09T35fM85N1WFJKktvzTuAiRJS89wl6QGGe6S1CDDXZIaZLhLUoNOGXcBAGeffXZNTk6OuwxJWlF2797946qa6LXuFRHuk5OTTE9Pj7sMSVpRkvzXidY5LCNJDTLcJalBhrskNchwl6QG9Q33JK9K8nCSR5PsS/Kxrv2zSX6QZE83berak+TTSQ4k2ZvkglGfhCTppQa5W+YF4JKqej7JqcADSf6lW/cXVfXF4/q/E9jYTW8Bbu0+JUnLpO+Ve815vls8tZtO9irJzcAd3XYPAmckWTt8qZKkQQ005p5kTZI9wBFgV1U91K26qRt6uSXJaV3bOuDJeZvPdG3H73NLkukk07Ozs0OcgiTpeAOFe1UdrapNwHrgwiRvAm4Efgf4XeAs4K+67um1ix773FZVU1U1NTHR8wErSdIiLegJ1ap6Jsk3gMuq6hNd8wtJ/gH48255Btgwb7P1wKFhCz2Rya33jGrXfR28+YqxHVuSTmaQu2UmkpzRzb8aeBvw3WPj6EkCXAk81m2yE3hvd9fMRcCzVXV4JNVLknoa5Mp9LbA9yRrmvgx2VNVXknw9yQRzwzB7gD/t+t8LXA4cAH4GvG/py5YknUzfcK+qvcD5PdovOUH/Aq4fvjRJ0mL5hKokNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSg/qGe5JXJXk4yaNJ9iX5WNd+XpKHknwvyReS/HLXflq3fKBbPznaU5AkHW+QK/cXgEuq6s3AJuCyJBcBHwduqaqNwE+A67r+1wE/qarfBG7p+kmSllHfcK85z3eLp3ZTAZcAX+zatwNXdvObu2W69ZcmyZJVLEnqa6Ax9yRrkuwBjgC7gP8EnqmqF7suM8C6bn4d8CRAt/5Z4Nd67HNLkukk07Ozs8OdhSTpJQYK96o6WlWbgPXAhcDre3XrPntdpdfLGqq2VdVUVU1NTEwMWq8kaQALulumqp4BvgFcBJyR5JRu1XrgUDc/A2wA6Nb/KvD0UhQrSRrMIHfLTCQ5o5t/NfA2YD9wP/BHXbdrgbu7+Z3dMt36r1fVy67cJUmjc0r/LqwFtidZw9yXwY6q+kqSx4E7k/w18G3gtq7/bcA/JjnA3BX71SOoW5J0En3Dvar2Auf3aP8+c+Pvx7f/D3DVklQnSVoUn1CVpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KD+oZ7kg1J7k+yP8m+JDd07R9N8qMke7rp8nnb3JjkQJInkrxjlCcgSXq5Uwbo8yLwoap6JMlrgd1JdnXrbqmqT8zvnOQNwNXAG4FfB/4tyW9V1dGlLFySdGJ9r9yr6nBVPdLNPwfsB9adZJPNwJ1V9UJV/QA4AFy4FMVKkgazoDH3JJPA+cBDXdMHkuxNcnuSM7u2dcCT8zab4eRfBpKkJTZwuCd5DXAX8MGq+ilwK/A6YBNwGPjksa49Nq8e+9uSZDrJ9Ozs7IILlySd2EDhnuRU5oL9c1X1JYCqeqqqjlbVz4HP8Iuhlxlgw7zN1wOHjt9nVW2rqqmqmpqYmBjmHCRJxxnkbpkAtwH7q+pT89rXzuv2buCxbn4ncHWS05KcB2wEHl66kiVJ/Qxyt8zFwHuA7yTZ07V9GLgmySbmhlwOAu8HqKp9SXYAjzN3p8313ikjScurb7hX1QP0Hke/9yTb3ATcNERdkqQh+ISqJDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoP6hnuSDUnuT7I/yb4kN3TtZyXZleR73eeZXXuSfDrJgSR7k1ww6pOQJL3UIFfuLwIfqqrXAxcB1yd5A7AVuK+qNgL3dcsA7wQ2dtMW4NYlr1qSdFJ9w72qDlfVI938c8B+YB2wGdjeddsOXNnNbwbuqDkPAmckWbvklUuSTmhBY+5JJoHzgYeAc6vqMMx9AQDndN3WAU/O22ymazt+X1uSTCeZnp2dXXjlkqQTGjjck7wGuAv4YFX99GRde7TVyxqqtlXVVFVNTUxMDFqGJGkAA4V7klOZC/bPVdWXuuanjg23dJ9HuvYZYMO8zdcDh5amXEnSIAa5WybAbcD+qvrUvFU7gWu7+WuBu+e1v7e7a+Yi4NljwzeSpOVxygB9LgbeA3wnyZ6u7cPAzcCOJNcBPwSu6tbdC1wOHAB+BrxvSSuWJPXVN9yr6gF6j6MDXNqjfwHXD1mXJGkIPqEqSQ0y3CWpQYa7JDVokB9UdQKTW+8Zy3EP3nzFWI4raeUw3FegcX2pgF8s0krhsIwkNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUF9wz3J7UmOJHlsXttHk/woyZ5uunzeuhuTHEjyRJJ3jKpwSdKJDXLl/lngsh7tt1TVpm66FyDJG4CrgTd22/xdkjVLVawkaTB9w72qvgk8PeD+NgN3VtULVfUD4ABw4RD1SZIWYZgx9w8k2dsN25zZta0DnpzXZ6Zre5kkW5JMJ5menZ0dogxJ0vEWG+63Aq8DNgGHgU927enRt3rtoKq2VdVUVU1NTEwssgxJUi+LCveqeqqqjlbVz4HP8Iuhlxlgw7yu64FDw5UoSVqoRYV7krXzFt8NHLuTZidwdZLTkpwHbAQeHq5ESdJCndKvQ5LPA28Fzk4yA3wEeGuSTcwNuRwE3g9QVfuS7AAeB14Erq+qo6MpXZJ0In3Dvaqu6dF820n63wTcNExRkqTh+ISqJDXIcJekBhnuktSgvmPu0nyTW+8Zy3EP3nzFWI4rrVReuUtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KD+oZ7ktuTHEny2Ly2s5LsSvK97vPMrj1JPp3kQJK9SS4YZfGSpN4GuXL/LHDZcW1bgfuqaiNwX7cM8E5gYzdtAW5dmjIlSQvRN9yr6pvA08c1bwa2d/PbgSvntd9Rcx4EzkiydqmKlSQNZrFj7udW1WGA7vOcrn0d8OS8fjNd28sk2ZJkOsn07OzsIsuQJPWy1D+opkdb9epYVduqaqqqpiYmJpa4DEla3RYb7k8dG27pPo907TPAhnn91gOHFl+eJGkxFhvuO4Fru/lrgbvntb+3u2vmIuDZY8M3kqTlc0q/Dkk+D7wVODvJDPAR4GZgR5LrgB8CV3Xd7wUuBw4APwPeN4KaJUl99A33qrrmBKsu7dG3gOuHLUqSNByfUJWkBhnuktQgw12SGmS4S1KD+v6gKr0STG69Z2zHPnjzFWM7trRYXrlLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkho01B/rSHIQeA44CrxYVVNJzgK+AEwCB4E/rqqfDFemJGkhluLK/feralNVTXXLW4H7qmojcF+3LElaRqMYltkMbO/mtwNXjuAYkqSTGDbcC/hakt1JtnRt51bVYYDu85xeGybZkmQ6yfTs7OyQZUiS5hv2D2RfXFWHkpwD7Ery3UE3rKptwDaAqampGrIOSdI8Q125V9Wh7vMI8GXgQuCpJGsBus8jwxYpSVqYRYd7ktOTvPbYPPB24DFgJ3Bt1+1a4O5hi5QkLcwwwzLnAl9Ocmw//1xVX03yLWBHkuuAHwJXDV+mJGkhFh3uVfV94M092v8buHSYoiRJwxn2B1WpeZNb7xnLcQ/efMVYjqs2+PoBSWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ3y9QPSK5SvPdAwvHKXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGjSzck1yW5IkkB5JsHdVxJEkvN5LXDyRZA/wt8AfADPCtJDur6vFRHE/S0hnXaw/AVx8spVG9W+ZC4EBVfR8gyZ3AZsBwl3RC4/xiGZdRfaGNKtzXAU/OW54B3jK/Q5ItwJZu8fkkTyzyWGcDP17ktiud5776rNbzhkbPPR/v2+Vk5/0bJ9poVOGeHm31koWqbcC2oQ+UTFfV1LD7WYk899V37qv1vGH1nvtiz3tUP6jOABvmLa8HDo3oWJKk44wq3L8FbExyXpJfBq4Gdo7oWJKk44xkWKaqXkzyAeBfgTXA7VW1bxTHYgmGdlYwz331Wa3nDav33Bd13qmq/r0kSSuKT6hKUoMMd0lq0IoO99X6ioMktyc5kuSxcdeynJJsSHJ/kv1J9iW5Ydw1LZckr0rycJJHu3P/2LhrWk5J1iT5dpKvjLuW5ZTkYJLvJNmTZHpB267UMffuFQf/wbxXHADXrIZXHCT5PeB54I6qetO461kuSdYCa6vqkSSvBXYDV66Sf/MAp1fV80lOBR4AbqiqB8dc2rJI8mfAFPArVfWucdezXJIcBKaqasEPb63kK/f/f8VBVf0vcOwVB82rqm8CT4+7juVWVYer6pFu/jlgP3NPQzev5jzfLZ7aTSvzymyBkqwHrgD+fty1rCQrOdx7veJgVfxHFySZBM4HHhpvJcunG5rYAxwBdlXVajn3vwH+Evj5uAsZgwK+lmR398qWga3kcO/7igO1KclrgLuAD1bVT8ddz3KpqqNVtYm5J74vTNL8kFySdwFHqmr3uGsZk4ur6gLgncD13ZDsQFZyuPuKg1WoG2++C/hcVX1p3PWMQ1U9A3wDuGzMpSyHi4E/7Mae7wQuSfJP4y1p+VTVoe7zCPBl5oajB7KSw91XHKwy3Y+KtwH7q+pT465nOSWZSHJGN/9q4G3Ad8db1ehV1Y1Vtb6qJpn7P/71qvqTMZe1LJKc3t04QJLTgbcDA98ht2LDvapeBI694mA/sGOErzh4RUnyeeDfgd9OMpPkunHXtEwuBt7D3NXbnm66fNxFLZO1wP1J9jJ3YbOrqlbVbYGr0LnAA0keBR4G7qmqrw668Yq9FVKSdGIr9spdknRihrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lq0P8Bvz499r3YUY4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "I = 1_000\n",
    "%time plt.hist(nonnegative_optimize(np.random.normal(1,1,size=I), I))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 7.48 s, sys: 160 ms, total: 7.64 s\n",
      "Wall time: 33.2 s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(array([3.5899e+04, 2.1123e+04, 1.9411e+04, 1.2947e+04, 7.0030e+03,\n",
       "        2.6140e+03, 7.7300e+02, 2.0100e+02, 2.4000e+01, 5.0000e+00]),\n",
       " array([0.        , 0.53820915, 1.07641831, 1.61462746, 2.15283661,\n",
       "        2.69104576, 3.22925492, 3.76746407, 4.30567322, 4.84388238,\n",
       "        5.38209153]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD4CAYAAAAO9oqkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAATXklEQVR4nO3df4xd5X3n8fenBhKUNGsTJsiyrTXqWm0cpBoyaywhrbIkawxUa1dKJJAarAjJbWSkRFu1Mf2H5gcS+aNhhZQguYsXs5uNY+WHsIJT1yJEEVIAD4ljMA71lHjD1BaerIGAohKZfvvHfby9Mnc81/Prjsfvl3R0z/2e55zzHCH4zHnOcy6pKiRJF7ffGXQHJEmDZxhIkgwDSZJhIEnCMJAkAZcMugNTdeWVV9bKlSsH3Q1JuqA8++yzv6qqobPrF2wYrFy5kpGRkUF3Q5IuKEn+b6+6w0SSJMNAkmQYSJIwDCRJGAaSJAwDSRKGgSQJw0CShGEgSeICfgN5OlZue2wg5z12360DOa8kTcY7A0nS5GGQ5N1JnknysySHk3y+1R9O8oskB9uyptWT5IEko0kOJbmu61ibkxxty+au+oeTPNf2eSBJZuNiJUm99TNM9BZwY1W9meRS4Mkk32/b/qKqvnVW+5uBVW25HngQuD7JFcA9wDBQwLNJ9lTVq63NFuApYC+wAfg+kqQ5MemdQXW82b5e2pY6xy4bgUfafk8Bi5MsBW4C9lfVqRYA+4ENbdv7qurHVVXAI8CmaVyTJOk89fXMIMmiJAeBk3T+g/5023RvGwq6P8m7Wm0Z8HLX7mOtdq76WI96r35sSTKSZGR8fLyfrkuS+tBXGFTV21W1BlgOrE1yDXA38AfAfwSuAD7Xmvca768p1Hv1Y3tVDVfV8NDQO/7fDJKkKTqv2URV9RrwQ2BDVZ1oQ0FvAf8TWNuajQErunZbDhyfpL68R12SNEf6mU00lGRxW78c+Bjw8zbWT5v5swl4vu2yB7ijzSpaB7xeVSeAfcD6JEuSLAHWA/vatjeSrGvHugN4dGYvU5J0Lv3MJloK7EyyiE547K6q7yX5QZIhOsM8B4E/a+33ArcAo8BvgE8BVNWpJF8EDrR2X6iqU23908DDwOV0ZhE5k0iS5tCkYVBVh4Bre9RvnKB9AVsn2LYD2NGjPgJcM1lfJEmzwzeQJUmGgSTJMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRKGgSQJw0CShGEgScIwkCRhGEiS6CMMkrw7yTNJfpbkcJLPt/rVSZ5OcjTJN5Nc1urvat9H2/aVXce6u9VfTHJTV31Dq40m2TbzlylJOpd+7gzeAm6sqj8E1gAbkqwDvgzcX1WrgFeBO1v7O4FXq+o/APe3diRZDdwGfAjYAHwtyaIki4CvAjcDq4HbW1tJ0hyZNAyq48329dK2FHAj8K1W3wlsausb23fa9o8mSavvqqq3quoXwCiwti2jVfVSVf0W2NXaSpLmSF/PDNpf8AeBk8B+4B+B16rqdGsyBixr68uAlwHa9teB93fXz9pnonqvfmxJMpJkZHx8vJ+uS5L60FcYVNXbVbUGWE7nL/kP9mrWPjPBtvOt9+rH9qoarqrhoaGhyTsuSerLec0mqqrXgB8C64DFSS5pm5YDx9v6GLACoG3/d8Cp7vpZ+0xUlyTNkX5mEw0lWdzWLwc+BhwBngA+3pptBh5t63vad9r2H1RVtfptbbbR1cAq4BngALCqzU66jM5D5j0zcXGSpP5cMnkTlgI726yf3wF2V9X3krwA7EryJeCnwEOt/UPA/0oySueO4DaAqjqcZDfwAnAa2FpVbwMkuQvYBywCdlTV4Rm7QknSpCYNg6o6BFzbo/4SnecHZ9f/GfjEBMe6F7i3R30vsLeP/kqSZoFvIEuSDANJkmEgScIwkCRhGEiSMAwkSRgGkiQMA0kShoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJoo8wSLIiyRNJjiQ5nOQzrf7XSf4pycG23NK1z91JRpO8mOSmrvqGVhtNsq2rfnWSp5McTfLNJJfN9IVKkibWz53BaeDPq+qDwDpga5LVbdv9VbWmLXsB2rbbgA8BG4CvJVmUZBHwVeBmYDVwe9dxvtyOtQp4Fbhzhq5PktSHScOgqk5U1U/a+hvAEWDZOXbZCOyqqreq6hfAKLC2LaNV9VJV/RbYBWxMEuBG4Ftt/53ApqlekCTp/J3XM4MkK4Frgadb6a4kh5LsSLKk1ZYBL3ftNtZqE9XfD7xWVafPqkuS5kjfYZDkvcC3gc9W1a+BB4HfA9YAJ4C/OdO0x+41hXqvPmxJMpJkZHx8vN+uS5Im0VcYJLmUThB8vaq+A1BVr1TV21X1L8Df0hkGgs5f9iu6dl8OHD9H/VfA4iSXnFV/h6raXlXDVTU8NDTUT9clSX3oZzZRgIeAI1X1la760q5mfww839b3ALcleVeSq4FVwDPAAWBVmzl0GZ2HzHuqqoAngI+3/TcDj07vsiRJ5+OSyZtwA/BJ4LkkB1vtr+jMBlpDZ0jnGPCnAFV1OMlu4AU6M5G2VtXbAEnuAvYBi4AdVXW4He9zwK4kXwJ+Sid8JElzZNIwqKon6T2uv/cc+9wL3NujvrfXflX1Ev82zCRJmmO+gSxJMgwkSYaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRKGgSQJw0CShGEgScIwkCRhGEiSMAwkSRgGkiT6+99eaoas3PbYwM597L5bB3ZuSfOfdwaSJMNAktRHGCRZkeSJJEeSHE7ymVa/Isn+JEfb55JWT5IHkowmOZTkuq5jbW7tjybZ3FX/cJLn2j4PJMlsXKwkqbd+7gxOA39eVR8E1gFbk6wGtgGPV9Uq4PH2HeBmYFVbtgAPQic8gHuA64G1wD1nAqS12dK134bpX5okqV+ThkFVnaiqn7T1N4AjwDJgI7CzNdsJbGrrG4FHquMpYHGSpcBNwP6qOlVVrwL7gQ1t2/uq6sdVVcAjXceSJM2B83pmkGQlcC3wNHBVVZ2ATmAAH2jNlgEvd+021mrnqo/1qPc6/5YkI0lGxsfHz6frkqRz6DsMkrwX+Dbw2ar69bma9qjVFOrvLFZtr6rhqhoeGhqarMuSpD71FQZJLqUTBF+vqu+08ittiIf2ebLVx4AVXbsvB45PUl/eoy5JmiP9zCYK8BBwpKq+0rVpD3BmRtBm4NGu+h1tVtE64PU2jLQPWJ9kSXtwvB7Y17a9kWRdO9cdXceSJM2Bft5AvgH4JPBckoOt9lfAfcDuJHcCvwQ+0bbtBW4BRoHfAJ8CqKpTSb4IHGjtvlBVp9r6p4GHgcuB77dFM2hQbz/75rN0YZg0DKrqSXqP6wN8tEf7ArZOcKwdwI4e9RHgmsn6IkmaHb6BLEkyDCRJhoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRKGgSQJw0CSRB9hkGRHkpNJnu+q/XWSf0pysC23dG27O8lokheT3NRV39Bqo0m2ddWvTvJ0kqNJvpnkspm8QEnS5Pq5M3gY2NCjfn9VrWnLXoAkq4HbgA+1fb6WZFGSRcBXgZuB1cDtrS3Al9uxVgGvAndO54IkSedv0jCoqh8Bp/o83kZgV1W9VVW/AEaBtW0ZraqXquq3wC5gY5IANwLfavvvBDad5zVIkqZpOs8M7kpyqA0jLWm1ZcDLXW3GWm2i+vuB16rq9Fn1npJsSTKSZGR8fHwaXZckdZtqGDwI/B6wBjgB/E2rp0fbmkK9p6raXlXDVTU8NDR0fj2WJE3okqnsVFWvnFlP8rfA99rXMWBFV9PlwPG23qv+K2Bxkkva3UF3e0nSHJnSnUGSpV1f/xg4M9NoD3BbkncluRpYBTwDHABWtZlDl9F5yLynqgp4Avh4238z8OhU+iRJmrpJ7wySfAP4CHBlkjHgHuAjSdbQGdI5BvwpQFUdTrIbeAE4DWytqrfbce4C9gGLgB1Vdbid4nPAriRfAn4KPDRjVydJ6sukYVBVt/coT/gf7Kq6F7i3R30vsLdH/SU6s40kSQPiG8iSJMNAkmQYSJIwDCRJGAaSJAwDSRJTfANZ6tfKbY8N7NzH7rt1YOeWLjTeGUiSDANJkmEgScIwkCRhGEiSMAwkSRgGkiQMA0kShoEkCcNAkoRhIEnCMJAkYRhIkugjDJLsSHIyyfNdtSuS7E9ytH0uafUkeSDJaJJDSa7r2mdza380yeau+oeTPNf2eSBJZvoiJUnn1s+dwcPAhrNq24DHq2oV8Hj7DnAzsKotW4AHoRMewD3A9cBa4J4zAdLabOna7+xzSZJm2aRhUFU/Ak6dVd4I7GzrO4FNXfVHquMpYHGSpcBNwP6qOlVVrwL7gQ1t2/uq6sdVVcAjXceSJM2RqT4zuKqqTgC0zw+0+jLg5a52Y612rvpYj3pPSbYkGUkyMj4+PsWuS5LONtMPkHuN99cU6j1V1faqGq6q4aGhoSl2UZJ0tqmGwSttiIf2ebLVx4AVXe2WA8cnqS/vUZckzaGphsEe4MyMoM3Ao131O9qsonXA620YaR+wPsmS9uB4PbCvbXsjybo2i+iOrmNJkubIJZM1SPIN4CPAlUnG6MwKug/YneRO4JfAJ1rzvcAtwCjwG+BTAFV1KskXgQOt3Req6sxD6U/TmbF0OfD9tkiS5tCkYVBVt0+w6aM92hawdYLj7AB29KiPANdM1g9J0uzxDWRJkmEgSTIMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJ9PHbRNKFauW2xwZy3mP33TqQ80rT4Z2BJMkwkCQZBpIkDANJEoaBJAnDQJKEYSBJwjCQJDHNMEhyLMlzSQ4mGWm1K5LsT3K0fS5p9SR5IMlokkNJrus6zubW/miSzdO7JEnS+ZqJO4P/XFVrqmq4fd8GPF5Vq4DH23eAm4FVbdkCPAid8ADuAa4H1gL3nAkQSdLcmI1hoo3Azra+E9jUVX+kOp4CFidZCtwE7K+qU1X1KrAf2DAL/ZIkTWC6YVDA3yd5NsmWVruqqk4AtM8PtPoy4OWufcdabaL6OyTZkmQkycj4+Pg0uy5JOmO6P1R3Q1UdT/IBYH+Sn5+jbXrU6hz1dxartgPbAYaHh3u2kSSdv2ndGVTV8fZ5EvgunTH/V9rwD+3zZGs+Bqzo2n05cPwcdUnSHJlyGCR5T5LfPbMOrAeeB/YAZ2YEbQYebet7gDvarKJ1wOttGGkfsD7JkvbgeH2rSZLmyHSGia4CvpvkzHH+T1X9XZIDwO4kdwK/BD7R2u8FbgFGgd8AnwKoqlNJvggcaO2+UFWnptEvSdJ5mnIYVNVLwB/2qP8/4KM96gVsneBYO4AdU+2LJGl6fANZkmQYSJIMA0kShoEkCcNAkoRhIEnCMJAkYRhIkpj+D9VJOsvKbY8N7NzH7rt1YOfWhc07A0mSYSBJMgwkSRgGkiQMA0kShoEkCcNAkoRhIEnCMJAk4RvI0oIyqLefffP5wuedgSRp/oRBkg1JXkwymmTboPsjSReTeREGSRYBXwVuBlYDtydZPdheSdLFY748M1gLjFbVSwBJdgEbgRcG2itJffFZxYVvvoTBMuDlru9jwPVnN0qyBdjSvr6Z5MUpnu9K4FdT3PdCsNCvD7zGhWDa15cvz1BPZs98/Gf473sV50sYpEet3lGo2g5sn/bJkpGqGp7ucearhX594DUuBAv9+uDCusZ58cyAzp3Aiq7vy4HjA+qLJF105ksYHABWJbk6yWXAbcCeAfdJki4a82KYqKpOJ7kL2AcsAnZU1eFZPOW0h5rmuYV+feA1LgQL/frgArrGVL1jaF6SdJGZL8NEkqQBMgwkSRdXGCz0n7xIsiPJySTPD7ovsyXJiiRPJDmS5HCSzwy6TzMpybuTPJPkZ+36Pj/oPs2GJIuS/DTJ9wbdl9mQ5FiS55IcTDIy6P7046J5ZtB+8uIfgP9CZyrrAeD2qlowbzkn+U/Am8AjVXXNoPszG5IsBZZW1U+S/C7wLLBpofxzTBLgPVX1ZpJLgSeBz1TVUwPu2oxK8t+AYeB9VfVHg+7PTEtyDBiuqvn2wtmELqY7g///kxdV9VvgzE9eLBhV9SPg1KD7MZuq6kRV/aStvwEcofMG+4JQHW+2r5e2ZUH9xZZkOXAr8D8G3Rf9m4spDHr95MWC+Y/IxSjJSuBa4OnB9mRmtSGUg8BJYH9VLajrA/478JfAvwy6I7OogL9P8mz7GZ1572IKg75+8kIXhiTvBb4NfLaqfj3o/sykqnq7qtbQeRN/bZIFM+SX5I+Ak1X17KD7MstuqKrr6PwS89Y2hDuvXUxh4E9eLBBtLP3bwNer6juD7s9sqarXgB8CGwbclZl0A/Bf25j6LuDGJP97sF2aeVV1vH2eBL5LZ5h6XruYwsCfvFgA2gPWh4AjVfWVQfdnpiUZSrK4rV8OfAz4+WB7NXOq6u6qWl5VK+n8O/iDqvqTAXdrRiV5T5vcQJL3AOuBeT/D76IJg6o6DZz5yYsjwO5Z/smLOZfkG8CPgd9PMpbkzkH3aRbcAHySzl+UB9tyy6A7NYOWAk8kOUTnD5j9VbUgp18uYFcBTyb5GfAM8FhV/d2A+zSpi2ZqqSRpYhfNnYEkaWKGgSTJMJAkGQaSJAwDSRKGgSQJw0CSBPwrSwoCWwKYEBsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "I = 100_000\n",
    "%time plt.hist(nonnegative_optimize(np.random.normal(1,1,size=I), I))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "dismod_mr",
   "language": "python",
   "name": "dismod_mr"
  },
  "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}