nonlinear_functions.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CP-SAT: Non-linear functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ortools.sat.python import cp_model\n",
    "\n",
    "import numpy as np\n",
    "from typing import Dict\n",
    "from math import exp, log\n",
    "\n",
    "FLOAT_APPROX_PRECISION = 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils import create_boolean_is_equal_to"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CP-SAT doesn't want you to use non-linear functions"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "CP SAT supports many integer operations. But you can't compose an IntVar by a custom function, in particular if it's non-linear. \n",
    "\n",
    "For example, you can't do exp(var)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "must be real number, not IntVar",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m/home/haxxor/projects/cpsat-fun/nonlinear_functions.ipynb Cell 6\u001b[0m in \u001b[0;36m<cell line: 8>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      <a href='vscode-notebook-cell:/home/haxxor/projects/cpsat-fun/nonlinear_functions.ipynb#W3sZmlsZQ%3D%3D?line=3'>4</a>\u001b[0m result \u001b[39m=\u001b[39m model\u001b[39m.\u001b[39mNewIntVar(\u001b[39m0\u001b[39m, \u001b[39m100\u001b[39m \u001b[39m*\u001b[39m \u001b[39m100\u001b[39m, \u001b[39m\"\u001b[39m\u001b[39mresult\u001b[39m\u001b[39m\"\u001b[39m)\n\u001b[1;32m      <a href='vscode-notebook-cell:/home/haxxor/projects/cpsat-fun/nonlinear_functions.ipynb#W3sZmlsZQ%3D%3D?line=5'>6</a>\u001b[0m \u001b[39m# Let's say we want result == exp(x)\u001b[39;00m\n\u001b[1;32m      <a href='vscode-notebook-cell:/home/haxxor/projects/cpsat-fun/nonlinear_functions.ipynb#W3sZmlsZQ%3D%3D?line=6'>7</a>\u001b[0m \u001b[39m# This will fail (rightfully)\u001b[39;00m\n\u001b[0;32m----> <a href='vscode-notebook-cell:/home/haxxor/projects/cpsat-fun/nonlinear_functions.ipynb#W3sZmlsZQ%3D%3D?line=7'>8</a>\u001b[0m model\u001b[39m.\u001b[39mAdd(result \u001b[39m==\u001b[39m exp(x))\n",
      "\u001b[0;31mTypeError\u001b[0m: must be real number, not IntVar"
     ]
    }
   ],
   "source": [
    "model = cp_model.CpModel()\n",
    "\n",
    "x = model.NewIntVar(0, 100, \"x\")\n",
    "result = model.NewIntVar(0, 100 * 100, \"result\")\n",
    "\n",
    "# Let's say we want result == exp(x)\n",
    "# This will fail (rightfully)\n",
    "model.Add(result == exp(x))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Precomputing the function values\n",
    "\n",
    "Since integer variables take a finite amount of values, we can simply pre-compute the non-linear function values, and store them in a map.\n",
    "\n",
    "Then, we set the image value to be equal to the value in the mapping when needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lookup_value_in_dict(\n",
    "    model: cp_model.CpModel,\n",
    "    key_var: cp_model.IntVar,\n",
    "    mapping: Dict[int, int],\n",
    "    mapping_name: str = \"mapping\",\n",
    "):\n",
    "    \"\"\"Creates a new variable equals to mapping[key_var]\n",
    "    If the value is absent from the mapping, the value_var takes the value 0.\n",
    "    \"\"\"\n",
    "    value_var = model.NewIntVar(\n",
    "        min(mapping.values()), max(mapping.values()), f\"{mapping_name}_{key_var.Name()}\"\n",
    "    )\n",
    "    for mapping_key, mapping_value in mapping.items():\n",
    "        key_var_is_equal_to = create_boolean_is_equal_to(model, key_var, mapping_key)\n",
    "        model.Add(value_var == mapping_value).OnlyEnforceIf(key_var_is_equal_to)\n",
    "\n",
    "    return value_var"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Example with the exponential and natural logarithm function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 412,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO : Make this work with variable bounds\n",
    "\n",
    "def exp_of_x(\n",
    "    model: cp_model.CpModel,\n",
    "    var: cp_model.IntVar,\n",
    "    float_precision_var=FLOAT_APPROX_PRECISION,\n",
    "    float_precision_image=FLOAT_APPROX_PRECISION,\n",
    "):\n",
    "    lb = var.Proto().domain[0]\n",
    "    ub = var.Proto().domain[1]\n",
    "    x_to_exp_x = {\n",
    "        x: round(exp(x / float_precision_var) * float_precision_image)\n",
    "        for x in range(lb, ub + 1)\n",
    "    }\n",
    "    exp_of_var = lookup_value_in_dict(model, var, x_to_exp_x, mapping_name=\"exp\")\n",
    "    return exp_of_var\n",
    "\n",
    "def log_of_x(\n",
    "    model: cp_model.CpModel,\n",
    "    var: cp_model.IntVar,\n",
    "    float_precision_var=FLOAT_APPROX_PRECISION,\n",
    "    float_precision_image=FLOAT_APPROX_PRECISION,\n",
    "):\n",
    "    # Log is only defined for x > 0\n",
    "    lb = max(var.Proto().domain[0], 1)\n",
    "    ub = max(var.Proto().domain[1], 1)\n",
    "    assert lb <= ub\n",
    "    x_to_log_x = {\n",
    "        x: round(log(x / float_precision_var) * float_precision_image)\n",
    "        for x in range(lb, ub + 1)\n",
    "    }\n",
    "    log_of_var = lookup_value_in_dict(model, var, x_to_log_x, mapping_name=\"log\")\n",
    "    return log_of_var"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's how it would be used in a model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 413,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Status = OPTIMAL\n",
      "Solver value of exp(0.69) is: 1.9937 (Real value: 1.9937155332430823)\n",
      "Solver value of log(0.69) is: -0.3711 (Real value: -0.37106368139083207)\n"
     ]
    }
   ],
   "source": [
    "float_precision_var = 100\n",
    "float_precision_image = 10_000\n",
    "value_of_x = 0.69\n",
    "\n",
    "model = cp_model.CpModel()\n",
    "\n",
    "x = model.NewIntVar(-10 * float_precision_var, 10 * float_precision_var, \"x\")\n",
    "exp_x_var = exp_of_x(\n",
    "    model,\n",
    "    x,\n",
    "    float_precision_var=float_precision_var,\n",
    "    float_precision_image=float_precision_image,\n",
    ")\n",
    "log_x_var = log_of_x(\n",
    "    model,\n",
    "    x,\n",
    "    float_precision_var=float_precision_var,\n",
    "    float_precision_image=float_precision_image,\n",
    ")\n",
    "model.Add(x == round(float_precision_var * value_of_x))\n",
    "\n",
    "solver = cp_model.CpSolver()\n",
    "status = solver.Solve(model)\n",
    "\n",
    "print(f\"Status = {solver.StatusName(status)}\")\n",
    "\n",
    "print(\n",
    "    f\"Solver value of exp({value_of_x}) is: {solver.Value(exp_x_var)/ float_precision_image} \"\n",
    "    + f\"(Real value: {exp(value_of_x)})\"\n",
    ")\n",
    "print(\n",
    "    f\"Solver value of log({value_of_x}) is: {solver.Value(log_x_var)/ float_precision_image} \"\n",
    "    + f\"(Real value: {log(value_of_x)})\"\n",
    ")\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Application: sum of logarithms\n",
    "\n",
    "### Approximate product of large numbers"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Multiplying large int values can become difficult as we reach the upper limit fast. For example, we can't store the product of x and y below, because the result _could_ be higher than `2**63 - 1`, which is the upper bound for ints in CP-SAT. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 414,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "__init__(): incompatible constructor arguments. The following argument types are supported:\n    1. ortools.util.python.sorted_interval_list.Domain(arg0: int, arg1: int)\n\nInvoked with: 1, 100000000000000000000",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m/home/haxxor/projects/cpsat-fun/nonlinear_functions.ipynb Cell 15\u001b[0m in \u001b[0;36m<cell line: 7>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      <a href='vscode-notebook-cell:/home/haxxor/projects/cpsat-fun/nonlinear_functions.ipynb#X16sZmlsZQ%3D%3D?line=3'>4</a>\u001b[0m y \u001b[39m=\u001b[39m model\u001b[39m.\u001b[39mNewIntVar(\u001b[39m10\u001b[39m\u001b[39m*\u001b[39m\u001b[39m*\u001b[39m\u001b[39m12\u001b[39m \u001b[39m-\u001b[39m \u001b[39m1000\u001b[39m, \u001b[39m10\u001b[39m\u001b[39m*\u001b[39m\u001b[39m*\u001b[39m\u001b[39m12\u001b[39m \u001b[39m+\u001b[39m \u001b[39m1000\u001b[39m, \u001b[39m\"\u001b[39m\u001b[39my\u001b[39m\u001b[39m\"\u001b[39m)\n\u001b[1;32m      <a href='vscode-notebook-cell:/home/haxxor/projects/cpsat-fun/nonlinear_functions.ipynb#X16sZmlsZQ%3D%3D?line=5'>6</a>\u001b[0m \u001b[39m# This will fail because the upper bound is above 2**63 - 1\u001b[39;00m\n\u001b[0;32m----> <a href='vscode-notebook-cell:/home/haxxor/projects/cpsat-fun/nonlinear_functions.ipynb#X16sZmlsZQ%3D%3D?line=6'>7</a>\u001b[0m x_times_y \u001b[39m=\u001b[39m model\u001b[39m.\u001b[39;49mNewIntVar(\u001b[39m1\u001b[39;49m, \u001b[39m10\u001b[39;49m\u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49m\u001b[39m10\u001b[39;49m \u001b[39m*\u001b[39;49m \u001b[39m10\u001b[39;49m\u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49m\u001b[39m10\u001b[39;49m, \u001b[39m\"\u001b[39;49m\u001b[39my\u001b[39;49m\u001b[39m\"\u001b[39;49m)\n",
      "File \u001b[0;32m~/ynstant/ynstaenv/lib/python3.9/site-packages/ortools/sat/python/cp_model.py:999\u001b[0m, in \u001b[0;36mCpModel.NewIntVar\u001b[0;34m(self, lb, ub, name)\u001b[0m\n\u001b[1;32m    983\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mNewIntVar\u001b[39m(\u001b[39mself\u001b[39m, lb, ub, name):\n\u001b[1;32m    984\u001b[0m     \u001b[39m\"\"\"Create an integer variable with domain [lb, ub].\u001b[39;00m\n\u001b[1;32m    985\u001b[0m \n\u001b[1;32m    986\u001b[0m \u001b[39mThe CP-SAT solver is limited to integer variables. If you have fractional\u001b[39;00m\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    996\u001b[0m \u001b[39m  a variable whose domain is [lb, ub].\u001b[39;00m\n\u001b[1;32m    997\u001b[0m \u001b[39m\"\"\"\u001b[39;00m\n\u001b[0;32m--> 999\u001b[0m     \u001b[39mreturn\u001b[39;00m IntVar(\u001b[39mself\u001b[39m\u001b[39m.\u001b[39m__model, Domain(lb, ub), name)\n",
      "\u001b[0;31mTypeError\u001b[0m: __init__(): incompatible constructor arguments. The following argument types are supported:\n    1. ortools.util.python.sorted_interval_list.Domain(arg0: int, arg1: int)\n\nInvoked with: 1, 100000000000000000000"
     ]
    }
   ],
   "source": [
    "model = cp_model.CpModel()\n",
    "\n",
    "x = model.NewIntVar(10**12 - 1000, 10**12 + 1000, \"x\")\n",
    "y = model.NewIntVar(10**12 - 1000, 10**12 + 1000, \"y\")\n",
    "\n",
    "# This will fail because the upper bound is above 2**63 - 1\n",
    "x_times_y = model.NewIntVar(1, 10**10 * 10**10, \"y\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using a lower upper bound for the variable equal to the multiplication of the two others will fail. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 415,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Status = MODEL_INVALID\n"
     ]
    }
   ],
   "source": [
    "model = cp_model.CpModel()\n",
    "\n",
    "x = model.NewIntVar(10**12 - 1000, 10**12 + 1000, \"x\")\n",
    "y = model.NewIntVar(10**12 - 1000, 10**12 + 1000, \"y\")\n",
    "\n",
    "# If we set the \n",
    "x_times_y = model.NewIntVar(1, 2**63 - 1, \"y\")\n",
    "model.AddMultiplicationEquality(x_times_y, [x, y])\n",
    "\n",
    "solver = cp_model.CpSolver()\n",
    "status = solver.Solve(model)\n",
    "\n",
    "print(f\"Status = {solver.StatusName(status)}\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A trick is to use the logarithm property : `log(A) + log(B) = log(A*B)`. This allows us to deal with much smaller number. Indeed, `log(2**63-1) ~= 43.66`.\n",
    "\n",
    "However, since we are dealing with decimal approximations, we must add a lot of decimals to keep estimations precise enough. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 416,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Status = OPTIMAL\n",
      "Solver value of x*y is: 9.999977681453931e+23 (Real value: 1000000000006999999999970)\n",
      "Relative error: -2.23186160683225e-06\n"
     ]
    }
   ],
   "source": [
    "model = cp_model.CpModel()\n",
    "\n",
    "x = model.NewIntVar(10**12 - 1000, 10**12 + 1000, \"x\")\n",
    "y = model.NewIntVar(10**12 - 1000, 10**12 + 1000, \"y\")\n",
    "\n",
    "big_x_value = 10**12 + 10\n",
    "big_y_value = 10**12 - 3\n",
    "log_float_precision = 100_000\n",
    "\n",
    "# Example values of x and y\n",
    "model.Add(x == big_x_value)\n",
    "model.Add(y == big_y_value)\n",
    "\n",
    "# Compute the log. We set the precision to be low, as we deal with big values\n",
    "log_x_var = log_of_x(model, x, float_precision_var=1, float_precision_image=log_float_precision)\n",
    "log_y_var = log_of_x(model, y, float_precision_var=1, float_precision_image=log_float_precision)\n",
    "\n",
    "# If we set the \n",
    "x_times_y = model.NewIntVar(1, 2*round(log(10**12+1000))*log_float_precision, \"y\")\n",
    "model.Add(x_times_y == log_x_var + log_y_var)\n",
    "\n",
    "solver = cp_model.CpSolver()\n",
    "status = solver.Solve(model)\n",
    "\n",
    "print(f\"Status = {solver.StatusName(status)}\")\n",
    "print(\n",
    "    f\"Solver value of x*y is: {exp(solver.Value(x_times_y) / log_float_precision)} \"\n",
    "    + f\"(Real value: {big_x_value * big_y_value})\"\n",
    ")\n",
    "print(f\"Relative error: {(exp(solver.Value(x_times_y) / log_float_precision) - (big_x_value * big_y_value) )/ (big_x_value * big_y_value)}\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Despite the relative error being small, it's still huge in absolute value. It's an approximation. \n",
    "\n",
    "Note as well, that the support of the functions must be small enough, since we need to store the logarithm image for every value of the interval. \n",
    "\n",
    "A simpler way could be to compute `x//BIG_NUMBER` and `y//BIG_NUMBER`, then set `z == x//BIG_NUMBER * y//BIG_NUMBER `, and interpret it as `z*BIG_NUMBER**2`."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Large product of many, many variables\n",
    "\n",
    "You may not always get a large product of two single variables. But you might get a large products of numerous smaller variables. \n",
    "\n",
    "With the logarithm, we get for free product of multiple variables. That's because  and $\\sum_k \\log(x_k) = \\log (\\prod_k x_k )$\n",
    "\n",
    "Once again, we are able to manipulate large quantities with an ok-level of approximation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 428,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Set x_0=28 \n",
      "Set x_1=74 \n",
      "Set x_2=74 \n",
      "Set x_3=38 \n",
      "Set x_4=17 \n",
      "Set x_5=66 \n",
      "Set x_6=42 \n",
      "Set x_7=71 \n",
      "Set x_8=25 \n",
      "Set x_9=86 \n",
      "... (100 variables in total)\n",
      "Status = OPTIMAL\n",
      "Solution found by the solver is: 2.1290700900253401e+158\n",
      "Real solution is: 212899174811596181578448858435179265963736896985091695902063357071093077381865583024923475729084353144425894505849610061214380949858418688000000000000000000000\n",
      "Relative error: 3.679765759901138e-05\n"
     ]
    }
   ],
   "source": [
    "from random import randint\n",
    "\n",
    "model = cp_model.CpModel()\n",
    "\n",
    "float_precision_log = 100_000\n",
    "\n",
    "n_variables = 100\n",
    "many_variables = [model.NewIntVar(1, 100, f\"{x}_{i}\") for i in range(n_variables)]\n",
    "many_log_variables = []\n",
    "\n",
    "log_result = model.NewIntVar(\n",
    "    1, 2 * float_precision_log * round(log(100)) * n_variables, \"result\"\n",
    ")\n",
    "\n",
    "# Set the var to be equal to random int between 1 and 10\n",
    "real_product = 1\n",
    "for i, var in enumerate(many_variables):\n",
    "    random_value = randint(1, 100)\n",
    "    if i < 10:\n",
    "        print(f\"Set {x}_{i}={random_value} \")\n",
    "    model.Add(var == random_value)\n",
    "\n",
    "    log_random_value = log_of_x(\n",
    "        model, var, float_precision_var=1, float_precision_image=float_precision_log\n",
    "    )\n",
    "    many_log_variables.append(log_random_value)\n",
    "\n",
    "    real_product *= random_value\n",
    "\n",
    "print(f\"... ({len(many_variables)} variables in total)\")\n",
    "\n",
    "model.Add(log_result == sum(many_log_variables))\n",
    "\n",
    "solver = cp_model.CpSolver()\n",
    "status = solver.Solve(model)\n",
    "\n",
    "print(f\"Status = {solver.StatusName(status)}\")\n",
    "print(\n",
    "    f\"Solution found by the solver is: {exp(solver.Value(log_result) / float_precision_log )}\"\n",
    ")\n",
    "print(f\"Real solution is: {real_product}\")\n",
    "print(\n",
    "    f\"Relative error: {(exp(solver.Value(log_result) / float_precision_log ) - real_product)/real_product}\"\n",
    ")\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exponentiation\n",
    "\n",
    "We know that $x^y = e^{y \\log (x)}$. Hence, we can manipulate decimal powers in the solver, such as square root."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 457,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Status = OPTIMAL\n",
      "Solution found by the solver is: 6.083012194367712\n",
      "Real solution is: 6.082762530298219\n"
     ]
    }
   ],
   "source": [
    "from math import sqrt\n",
    "\n",
    "model = cp_model.CpModel()\n",
    "\n",
    "float_precision_log = 10_000\n",
    "value_of_x = 37\n",
    "\n",
    "x = model.NewIntVar(1, 100, \"x\")\n",
    "model.Add(x == value_of_x)\n",
    "\n",
    "log_x = log_of_x(\n",
    "    model, x, float_precision_var=1, float_precision_image=float_precision_log\n",
    ")\n",
    "\n",
    "sqrt_x = model.NewIntVar(1, 10 * float_precision_log * float_precision_log, \"sqrt_x\")\n",
    "# we want sqrt_x == 0.5 * log_x\n",
    "model.Add(float_precision_log * sqrt_x >= round(0.5 * float_precision_log) * log_x)\n",
    "\n",
    "solver = cp_model.CpSolver()\n",
    "status = solver.Solve(model)\n",
    "\n",
    "print(f\"Status = {solver.StatusName(status)}\")\n",
    "print(\n",
    "    f\"Solution found by the solver is: {exp(solver.Value(sqrt_x) / float_precision_log )}\"\n",
    ")\n",
    "print(f\"Real solution is: {sqrt(value_of_x)}\")\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ynstaenv",
   "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.9.5"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "60aad423cf61d8722f5f696b996614860486ab31702984425ac69059c5f0f73d"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}