{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e9906f9b-affc-4b8f-a9ea-3c7ab4045332",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import optimize\n",
    "import pdfo\n",
    "import x3cflux\n",
    "import hopsy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b80366c1-c056-4fab-a7c9-38d8a976f2b0",
   "metadata": {},
   "source": [
    "# Interoperation with Optimization Packages\n",
    "\n",
    "Via Python, a huge variety of optimization packages are available. Still, some restrictions have to be met as these software packages are usually implemented from a methodological point of view. *13CFLUX*, however, was designed to match this perspective and therefore seamlessly operates with many packages. This notebook demonstrates, how *13CFLUX* integrates into two optimization packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4eaa4f49-95c5-479f-80c9-f0bee66788f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "simulator = x3cflux.create_simulator_from_fml(\"spiralus.fml\", \"ms+uptake\")\n",
    "mle = x3cflux.get_parameters(simulator.parameter_space, simulator.configurations[0].parameter_entries)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "49f54ad3-c466-4f60-a73b-2eec69a19c22",
   "metadata": {},
   "outputs": [],
   "source": [
    "ineq_sys = simulator.parameter_space.inequality_system"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8f15a73-b6ca-4e80-83ce-9da8018cefce",
   "metadata": {},
   "source": [
    "Generate starting point and constraints using *hopsy*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b31b7608-fa32-425a-978e-44305b5d0536",
   "metadata": {},
   "outputs": [],
   "source": [
    "problem = hopsy.Problem(ineq_sys.matrix, ineq_sys.bound)\n",
    "problem = hopsy.add_box_constraints(problem, 2*[0], 2*[5])\n",
    "starting_point = hopsy.compute_chebyshev_center(problem)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b89f5a38-1d9c-441c-b192-d3864e8c12df",
   "metadata": {},
   "source": [
    "### Using solvers from *scipy*\n",
    "\n",
    "*scipy.minimize* is one of the most frequently used framework for general purpose optimization. With the interface functions and definition of linear inequality constraints, *13CFLUX3* integrates seamlessly. Here, *trust-constr* was used, but others are available to handle linear inequalities too (e.g. *SLSQP* and *COBYQA*)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "23347aa3-e634-41eb-9ee2-dd55b4f4087f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "| niter |f evals|CG iter|  obj func   |tr radius |   opt    |  c viol  |\n",
      "|-------|-------|-------|-------------|----------|----------|----------|\n",
      "|   1   |   1   |   0   | +3.2424e+03 | 1.00e+00 | 8.25e+02 | 0.00e+00 |\n",
      "|   2   |   2   |   1   | +1.8018e+03 | 2.00e+00 | 6.82e+02 | 0.00e+00 |\n",
      "|   3   |   3   |   3   | +9.4316e+02 | 9.23e+00 | 5.03e+02 | 0.00e+00 |\n",
      "|   4   |   4   |   5   | +6.2678e+02 | 9.23e+00 | 3.52e+02 | 0.00e+00 |\n",
      "|   5   |   5   |   7   | +4.4453e+02 | 9.23e+00 | 2.58e+02 | 0.00e+00 |\n",
      "|   6   |   6   |   9   | +3.1440e+01 | 1.01e+01 | 1.05e+02 | 0.00e+00 |\n",
      "|   7   |   7   |  11   | +3.1440e+01 | 1.01e+00 | 1.05e+02 | 0.00e+00 |\n",
      "|   8   |   8   |  13   | +1.8182e+01 | 1.01e+00 | 7.21e+01 | 0.00e+00 |\n",
      "|   9   |   9   |  15   | +1.8182e+01 | 2.85e-01 | 7.21e+01 | 0.00e+00 |\n",
      "|  10   |  10   |  17   | +3.1015e+00 | 1.23e+00 | 2.85e+01 | 0.00e+00 |\n",
      "|  11   |  11   |  19   | +5.2174e-01 | 1.23e+00 | 1.05e+01 | 0.00e+00 |\n",
      "|  12   |  12   |  21   | +1.2135e-02 | 1.23e+00 | 1.47e+00 | 0.00e+00 |\n",
      "|  13   |  13   |  23   | +1.3566e-03 | 1.23e+00 | 9.39e-02 | 0.00e+00 |\n",
      "|  14   |  14   |  25   | +4.2115e-04 | 1.23e+00 | 6.06e-03 | 0.00e+00 |\n",
      "|  15   |  14   |  25   | +4.2115e-04 | 6.13e+00 | 5.58e-02 | 0.00e+00 |\n",
      "|  16   |  15   |  27   | +1.6714e-05 | 6.13e+00 | 4.17e-03 | 0.00e+00 |\n",
      "|  17   |  15   |  27   | +1.6714e-05 | 3.06e+01 | 1.34e-02 | 0.00e+00 |\n",
      "|  18   |  16   |  29   | +6.3684e-07 | 3.06e+01 | 1.34e-04 | 0.00e+00 |\n",
      "|  19   |  16   |  29   | +6.3684e-07 | 1.53e+02 | 2.10e-03 | 0.00e+00 |\n",
      "|  20   |  17   |  31   | +2.5505e-08 | 1.53e+02 | 1.66e-05 | 0.00e+00 |\n",
      "|  21   |  17   |  31   | +2.5505e-08 | 7.66e+02 | 3.85e-04 | 0.00e+00 |\n",
      "|  22   |  18   |  33   | +1.0195e-09 | 7.66e+02 | 1.27e-06 | 0.00e+00 |\n",
      "\n",
      "`gtol` termination condition is satisfied.\n",
      "Number of iterations: 22, function evaluations: 18, CG iterations: 33, optimality: 1.27e-06, constraint violation: 0.00e+00, execution time: 0.13 s.\n"
     ]
    }
   ],
   "source": [
    "simulator.parameter_space.constraint_violation_tolerance = 1e-3\n",
    "x3cflux.omp.num_max_threads = 10\n",
    "result = optimize.minimize(simulator.compute_loss, starting_point.flatten(), jac=simulator.compute_loss_gradient,\n",
    "                           constraints=[optimize.LinearConstraint(ineq_sys.matrix, ub=ineq_sys.bound, keep_feasible=True)],\n",
    "                           method=\"trust-constr\", tol=1e-5, options={\"maxiter\": 50, \"verbose\": 2})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b3097c65-9bb6-4a48-b3ec-5dcbee04fb9c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.019462704647814e-09, 1.7322214344899365e-06)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result.fun, np.linalg.norm(result.x - mle)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb8b022f-60bd-4737-ada4-1018945933a6",
   "metadata": {},
   "source": [
    "### Using gradient-free solver *LINCOA*\n",
    "\n",
    "Gradient-based optimization is not always viable if gradient costs are high. This is especially the case for high numbers of parameters in large network models. *LINCOA* is one of Powel's derivative-free optimizers that is specifically designed for linear inequality constraints and is available via the Python package *pdfo*. Similar to the *scipy* example, *13CFLUX* integrates without notable overhead. Warnings are issued when the solver hits numerical edge-cases (here: evaluating at the boundary of the polytope), which are handeled automatically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b34d71a0-0d32-4c77-bf25-6b5e739b7ff5",
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[33mOperation \"SparseLES Solver\" failed: Singular system\u001b[0m\n",
      "\u001b[1m\u001b[33mOperation \"SparseLES Solver\" failed: Singular system\u001b[0m\n",
      "\u001b[1m\u001b[33mOperation \"SparseLES Solver\" failed: Singular system\u001b[0m\n",
      "\u001b[1m\u001b[33mOperation \"SparseLES Solver\" failed: Singular system\u001b[0m\n",
      "\u001b[1m\u001b[33mOperation \"SparseLES Solver\" failed: Singular system\u001b[0m\n",
      "\u001b[1m\u001b[33mOperation \"SparseLES Solver\" failed: Singular system\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "simulator.parameter_space.constraint_violation_tolerance = 1\n",
    "result = pdfo.pdfo(lambda x: simulator.compute_loss(x) if simulator.parameter_space.contains(x) else 1e20, \n",
    "                   starting_point.flatten(), \n",
    "                   constraints=optimize.LinearConstraint(ineq_sys.matrix, ub=ineq_sys.bound, lb=-np.infty),\n",
    "                   method=\"lincoa\", options={'maxfev': 20_000, 'radius_final': 1e-15})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "7d4a90eb",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2.0830858278492343e-28, 7.447602459741819e-16)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result.fun, np.linalg.norm(result.x - mle)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7ed22ca-0091-4174-b97b-69bc8a7bf635",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}