{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "d0996120-bb17-4476-b912-ce155100b2cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import minimize\n",
    "import numpy as np\n",
    "\n",
    "# Define Parameters\n",
    "lambda_vals = np.array([0.03, 0.04, 0.05, 0.06, 0.07, 1, 1.1, 1.2, 1.3, 1.4, 1.5])  # Request rates ascendingly\n",
    "N = len(lambda_vals)\n",
    "B = 4.4  # Cache size\n",
    "c_delta = 1  # Age linear cost\n",
    "c_f = 7  # Fetching linear cost (caching miss cost)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "fcf0c13c-5b2c-457e-9aa6-8d349fcf13fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def theoretical_opt(lambda_vals, B, c_f, c_delta):\n",
    "    \"\"\"\n",
    "    Perform theoretical optimization to compute optimal hit probabilities.\n",
    "    \n",
    "    Parameters:\n",
    "    - lambda_vals: Array of request rates.\n",
    "    - B: Cache size (constraint on total hit probabilities).\n",
    "    - c_f: Cost of fetching (cache miss cost).\n",
    "    - c_delta: Cost of caching (hit cost).\n",
    "\n",
    "    Returns:\n",
    "    - h_optt: Optimal hit probabilities.\n",
    "    \"\"\"\n",
    "    N = len(lambda_vals)\n",
    "    h_optt = np.zeros(N)  # Initialize optimal hit probabilities\n",
    "    differenc_set = np.arange(N)  # Set of variables to optimize\n",
    "    fix_i = []  # Set of fixed variables (those already optimized)\n",
    "    n = N\n",
    "    b = B\n",
    "    flag = True\n",
    "\n",
    "    while flag:\n",
    "        if n == 0:  # If no variables left to optimize\n",
    "            if b > 0:  # If there is leftover cache size, redistribute it\n",
    "                differenc_set = np.where(h_optt == 0)[0]  # Find zero hit probability variables\n",
    "                fix_i = np.setdiff1d(np.arange(N), differenc_set)\n",
    "                n = len(differenc_set)\n",
    "                continue\n",
    "            else:  # No variables to optimize, finalize\n",
    "                h_optt[differenc_set] = 0\n",
    "                break\n",
    "        \n",
    "        # Calculate the optimal Lagrange multiplier (mu) and hit probabilities for the set of variables\n",
    "        mu = max(0, (n * c_f - b * c_delta) / np.sum(1 / lambda_vals[differenc_set]))\n",
    "        h_optt[differenc_set] = (c_f - mu / lambda_vals[differenc_set]) / c_delta\n",
    "        \n",
    "        # If mu < 0, adjust the cache size to set mu to zero in the next iteration\n",
    "        if mu < 0:\n",
    "            b = (n * c_f / c_delta)\n",
    "            continue\n",
    "        \n",
    "        # Identify violations of the hit probability constraints (h > 1 or h < 0)\n",
    "        larger_i = np.where(h_optt > 1)[0]\n",
    "        smaller_i = np.where(h_optt < 0)[0]\n",
    "\n",
    "        # If no violations, the optimal solution is reached\n",
    "        if len(smaller_i) == 0 and len(larger_i) == 0:\n",
    "            break\n",
    "        \n",
    "        # Find the furthest object from the boundary (either 0 or 1)\n",
    "        min_viol = 0\n",
    "        min_viol_i = -1\n",
    "        if len(smaller_i) > 0:\n",
    "            min_viol, min_viol_i = np.min(h_optt[smaller_i]), np.argmin(h_optt[smaller_i])\n",
    "\n",
    "        max_viol = 0\n",
    "        max_viol_i = -1\n",
    "        if len(larger_i) > 0:\n",
    "            max_viol, max_viol_i = np.max(h_optt[larger_i] - 1), np.argmax(h_optt[larger_i] - 1)\n",
    "        \n",
    "        # Choose the variable with the largest violation to adjust\n",
    "        if max_viol > abs(min_viol):\n",
    "            viol_i = max_viol_i\n",
    "            min_viol_flag = 0\n",
    "        else:\n",
    "            viol_i = min_viol_i\n",
    "            min_viol_flag = 1\n",
    "        \n",
    "        # Set the furthest object to the nearest boundary (0 or 1)\n",
    "        if min_viol_flag:\n",
    "            h_optt[viol_i] = 0\n",
    "        else:\n",
    "            h_optt[viol_i] = min(1, b)\n",
    "        \n",
    "        # Update cache size and fix the selected variable\n",
    "        b -= h_optt[viol_i]\n",
    "        fix_i.append(viol_i)\n",
    "        differenc_set = np.setdiff1d(np.arange(N), fix_i)\n",
    "        n = N - len(fix_i)\n",
    "    \n",
    "    return h_optt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5d223324-d193-416a-b2c3-1e143e981a37",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def numerical_opt(lambda_vals, B, c_f, c_delta):\n",
    "    \"\"\"\n",
    "    Perform numerical optimization to compute optimal hit probabilities.\n",
    "\n",
    "    Parameters:\n",
    "    - lambda_vals: Array of request rates.\n",
    "    - B: Cache size (constraint on total hit probabilities).\n",
    "    - c_f: Cost of fetching (cache miss cost).\n",
    "    - c_delta: Cost of caching (hit cost).\n",
    "\n",
    "    Returns:\n",
    "    - x_opt: Optimal hit probabilities.\n",
    "    \"\"\"\n",
    "    N = len(lambda_vals)  # Number of items\n",
    "\n",
    "    # Initial guess: Even distribution of cache capacity\n",
    "    x_init = np.full(N, B / N)\n",
    "\n",
    "    # Objective function\n",
    "    def objective(x):\n",
    "        return np.sum(lambda_vals * ((1 - x) * c_f + x**2 * c_delta / 2))\n",
    "\n",
    "    # Constraint: Sum of hit probabilities <= cache size (B)\n",
    "    def constraint_total_hit(x):\n",
    "        return B - np.sum(x)  # Non-negative means constraint satisfied\n",
    "\n",
    "    # Bounds for hit probabilities: 0 <= h_i <= 1\n",
    "    bounds = [(0, 1) for _ in range(N)]\n",
    "\n",
    "    # Optimization\n",
    "    constraints = [{'type': 'ineq', 'fun': constraint_total_hit}]  # Inequality constraint\n",
    "    result = minimize(\n",
    "        objective, \n",
    "        x_init, \n",
    "        method='SLSQP',  # Sequential Least Squares Quadratic Programming\n",
    "        bounds=bounds, \n",
    "        constraints=constraints, \n",
    "        options={'disp': True}  # Set to True for optimization output\n",
    "    )\n",
    "\n",
    "    # Optimal solution\n",
    "    if result.success:\n",
    "        return result.x  # Optimal hit probabilities\n",
    "    else:\n",
    "        raise ValueError(\"Optimization failed: \" + result.message)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "073be740-dc97-454b-87e7-2f8f93c8f137",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully    (Exit mode 0)\n",
      "            Current function value: 16.15721739129548\n",
      "            Iterations: 4\n",
      "            Function evaluations: 48\n",
      "            Gradient evaluations: 4\n",
      "Hit probabilities (numerical and theoretical):\n",
      " [[0.00000000e+00 0.00000000e+00]\n",
      " [0.00000000e+00 0.00000000e+00]\n",
      " [0.00000000e+00 0.00000000e+00]\n",
      " [0.00000000e+00 0.00000000e+00]\n",
      " [0.00000000e+00 0.00000000e+00]\n",
      " [4.58923826e-13 1.00000000e+00]\n",
      " [4.26087031e-01 0.00000000e+00]\n",
      " [9.73912969e-01 4.00000000e-01]\n",
      " [1.00000000e+00 0.00000000e+00]\n",
      " [1.00000000e+00 0.00000000e+00]\n",
      " [1.00000000e+00 0.00000000e+00]]\n",
      "Objective function values (numerical, theoretical): [16.15721739129548, 44.486]\n",
      "Constraint violations (numerical, theoretical): [1.241673430740775e-12, -3.0]\n"
     ]
    }
   ],
   "source": [
    "\n",
    "# Optimization\n",
    "h_numerical = numerical_opt(lambda_vals, B, c_f, c_delta)\n",
    "h_theoretical = theoretical_opt(lambda_vals, B, c_f, c_delta)\n",
    "\n",
    "# Comparison of Hit Probabilities\n",
    "hit_opt = np.vstack((h_numerical, h_theoretical)).T  # Combine numerical and theoretical hit probabilities\n",
    "\n",
    "# Objective Function Calculation\n",
    "obj_1 = np.sum(lambda_vals * ((1 - h_numerical) * c_f + h_numerical**2 * c_delta / 2))\n",
    "obj_2 = np.sum(lambda_vals * ((1 - h_theoretical) * c_f + h_theoretical**2 * c_delta / 2))\n",
    "obj = [obj_1, obj_2]  # Store objective function values for both methods\n",
    "\n",
    "# Constraints\n",
    "const_1 = np.sum(h_numerical) - B\n",
    "const_2 = np.sum(h_theoretical) - B\n",
    "constraint = [const_1, const_2]  # Check if the cache size constraint is satisfied\n",
    "\n",
    "# Outputs\n",
    "print(\"Hit probabilities (numerical and theoretical):\\n\", hit_opt)\n",
    "print(\"Objective function values (numerical, theoretical):\", obj)\n",
    "print(\"Constraint violations (numerical, theoretical):\", constraint)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "graphs",
   "language": "python",
   "name": "graphs"
  },
  "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.12.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}