{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3eb0337f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ff6f670",
   "metadata": {},
   "outputs": [],
   "source": [
    "class InputError(Exception):\n",
    "    pass\n",
    "\n",
    "file_path = 'sol_sampling_lp.txt'\n",
    "def read_input(file_path):\n",
    "    try:\n",
    "        with open(file_path, 'r') as file:\n",
    "            lines = file.readlines()\n",
    "            if len(lines) < 12:\n",
    "                raise InputError(\"File has wrong format, expected at least 12 lines.\")\n",
    "            \n",
    "            if (lines[0].strip() != \"N\" or lines[2].strip() != \"sigma\" or lines[4].strip() != \"alpha\" or lines[6].strip() != \"beta\"):\n",
    "                raise InputError(\"File has wrong format, expected parameters N, sigma, alpha, beta.\")\n",
    "            split_line = lines[1].strip().split()\n",
    "            if (len(split_line) != 1):\n",
    "                raise InputError(\"File has wrong format, the second line should contain exactly the value of N.\")\n",
    "            N = int(lines[1].strip())\n",
    "            split_line = lines[3].strip().split()\n",
    "            if (len(split_line) != 1):\n",
    "                raise InputError(\"File has wrong format, the fourth line should contain exactly the value of sigma.\")\n",
    "            sigma = str(lines[3].strip())\n",
    "            split_line = lines[5].strip().split()\n",
    "            if (len(split_line) != 1):\n",
    "                raise InputError(\"File has wrong format, the sixth line should contain exactly the value of alpha.\")\n",
    "            alpha = str(lines[5].strip())\n",
    "            split_line = lines[7].strip().split()\n",
    "            if (len(split_line) != 1):\n",
    "                raise InputError(\"File has wrong format, the eighth line should contain exactly the value of beta.\")\n",
    "            beta = str(lines[7].strip())\n",
    "            \n",
    "            if (lines[8].strip() != \"y\"):\n",
    "                raise InputError(\"File has wrong format, the ninth line should be 'y'.\")\n",
    "            split_line = lines[9].strip().split()\n",
    "            if (len(split_line) != 1):\n",
    "                raise InputError(\"File has wrong format, the tenth line should contain exactly the value of y\")\n",
    "            y = int(split_line[0])\n",
    "            if (lines[10].strip() != \"x\"):\n",
    "                raise InputError(\"File has wrong format, the eleventh line should be 'x'.\")\n",
    "            i = 11\n",
    "            x_variables = {}\n",
    "            while (i < len(lines) and lines[i].strip() != \"v\"):\n",
    "                split_line = lines[i].strip().split()\n",
    "                if (len(split_line) != 3):\n",
    "                    raise InputError(\"File has wrong format, lines after 'x' but before 'v' should contain two indices and one value.\")\n",
    "                if int(split_line[0]) not in x_variables:\n",
    "                    x_variables[int(split_line[0])] = {}\n",
    "                x_variables[int(split_line[0])][int(split_line[1])] = int(split_line[2])\n",
    "                i += 1\n",
    "            # now lines[i] == \"v\"\n",
    "            i += 1\n",
    "            v_variables = {}\n",
    "            while (i < len(lines) and lines[i].strip() != \"w\"):\n",
    "                split_line = lines[i].strip().split()\n",
    "                if (len(split_line) != 3):\n",
    "                    raise InputError(\"File has wrong format, lines after 'v' should contain two indices and one value.\")\n",
    "                if int(split_line[0]) not in v_variables:\n",
    "                    v_variables[int(split_line[0])] = {}\n",
    "                v_variables[int(split_line[0])][int(split_line[1])] = int(split_line[2])\n",
    "                i += 1\n",
    "            # now lines[i] == \"w\"\n",
    "            i += 1\n",
    "            w_variables = {}\n",
    "            while (i < len(lines)):\n",
    "                split_line = lines[i].strip().split()\n",
    "                if (len(split_line) != 3):\n",
    "                    raise InputError(\"File has wrong format, lines after 'w' should contain two indices and one value.\")\n",
    "                if int(split_line[0]) not in w_variables:\n",
    "                    w_variables[int(split_line[0])] = {}\n",
    "                w_variables[int(split_line[0])][int(split_line[1])] = int(split_line[2])\n",
    "                i += 1\n",
    "            \n",
    "            print(\"Read input variables and parameters; parameters are: N = \" + str(N) + \", sigma = \" + sigma + \", alpha = \" + alpha + \", beta = \" + beta + \".\")\n",
    "            known_parameters = [[2500, \"0.663\", \"1.5\", \"0.01\"], [2500, \"0.623\", \"1\", \"0.01\"]]\n",
    "            if (not [N, sigma, alpha, beta] in known_parameters):\n",
    "                print (\"Warning: your input contains unexpected parameters for which we didn't manually precompute upper and lower bounds on certain expressions. We will compute these bounds on the fly, beware of rounding errors due to inexact floating point operations.\")\n",
    "            \n",
    "            return [y, x_variables, v_variables, w_variables], N, sigma, alpha, beta\n",
    "    except FileNotFoundError:\n",
    "        print(\"File not found.\")\n",
    "    except IOError:\n",
    "        print(\"An error occurred while reading the file.\")\n",
    "    except InputError as error:\n",
    "        print(\"InputError: \" + str(error))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f698159",
   "metadata": {},
   "outputs": [],
   "source": [
    "def value(variable_dict, *indices):\n",
    "    if indices[0] not in variable_dict:\n",
    "        return 0\n",
    "    if len(indices) == 1:\n",
    "        return variable_dict[indices[0]]\n",
    "    if indices[1] not in variable_dict[indices[0]]:\n",
    "        return 0\n",
    "    return variable_dict[indices[0]][indices[1]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "831dd1f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def upper_bound_e_to_minus_sigma_times_beta(sigma, beta): #scaled by 10^20\n",
    "    if (sigma == \"0.663\" and beta == \"0.01\"):\n",
    "        return 99339192995802757404\n",
    "    elif (sigma == \"0.623\" and beta == \"0.01\"):\n",
    "        return 99378936621196242200\n",
    "    else:\n",
    "        return int(math.exp(- float(sigma) * float(beta)) * pow(10, 20)) + 1\n",
    "def lower_bound_e_to_minus_one_half_beta(beta): #scaled by 10^20\n",
    "    if (beta == \"0.01\"):\n",
    "        return 99501247919268231335\n",
    "    else:\n",
    "        return int(math.exp(- float(beta)/2) * pow(10, 20))\n",
    "def lower_bound_e_to_minus_beta(beta): #scaled by 10^20\n",
    "    if (beta == \"0.01\"):\n",
    "        return 99004983374916805357\n",
    "    else:\n",
    "        return int(math.exp(- float(beta)) * pow(10, 20))\n",
    "def upper_bound_4_beta_e_to_sigma_beta(sigma, beta): #scaled by 10^20\n",
    "    if (sigma == \"0.663\" and beta == \"0.01\"):\n",
    "        return 4026608108411960334\n",
    "    elif (sigma == \"0.623\" and beta == \"0.01\"):\n",
    "        return 4024997787254298031\n",
    "    else:\n",
    "        return int(4 * float(beta) * math.exp(float(sigma) * float(beta)) * pow(10, 20)) + 1\n",
    "def upper_bound_delta_2(N, sigma, alpha, beta): #scaled by 10^20\n",
    "    if (N == 2500 and sigma == \"0.663\" and alpha == \"1.5\" and beta == \"0.01\"):\n",
    "        return 3811092096437300449\n",
    "    elif (N == 2500 and sigma == \"0.623\" and alpha == \"1\" and beta == \"0.01\"):\n",
    "        return 7375219460756672439\n",
    "    else:\n",
    "        factor1 = float(alpha) + 2 * float(beta) / (math.exp(N * float(sigma) * float(beta)) * (math.exp(float(sigma) * float(beta)) - 1))\n",
    "        factor2 = math.exp(-N * float(sigma) * float(beta)) / pow(1 - math.exp(-float(sigma) * float(beta)), 2)\n",
    "        factor3 = 1 + N - math.exp(-float(sigma) * float(beta)) * N\n",
    "        return int(factor1 * factor2 * factor3 * pow(10,20)) + 1\n",
    "def upper_bound_alpha_plus_delta_1_times_e_to_sigma_beta(N, sigma, alpha, beta):\n",
    "    if (N == 2500 and sigma == \"0.663\" and alpha == \"1.5\" and beta == \"0.01\"): #scaled by 10^20\n",
    "        return 150997842396816263968\n",
    "    elif (N == 2500 and sigma == \"0.623\" and alpha == \"1\" and beta == \"0.01\"):\n",
    "        return 100625055544583975096\n",
    "    else:\n",
    "        delta_1 = 4 * float(beta) / (math.exp(N * float(sigma) * float(beta)) * (math.exp(float(sigma) * float(beta)) - 1))\n",
    "        return int(((float(alpha) + delta_1) * math.exp(float(sigma) * float(beta))) * pow(10,20)) + 1\n",
    "def lower_bound_sigma_to_minus_2(sigma): #scaled by 10^10\n",
    "    if (sigma == \"0.663\"):\n",
    "        return 22749556952\n",
    "    elif (sigma == \"0.623\"):\n",
    "        return 25764629800\n",
    "    else:\n",
    "        return int(1 / (float(sigma) * float(sigma)) * pow(10,10))\n",
    "def upper_bound_sigma_squared(sigma): #scaled by 10^6\n",
    "    if (sigma == \"0.663\"):\n",
    "        return 439569\n",
    "    elif (sigma == \"0.623\"):\n",
    "        return 388129\n",
    "    else:\n",
    "        return int(float(sigma) * float(sigma) * pow(10,6)) + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b05362bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "class InfeasibilityError(Exception):\n",
    "    pass\n",
    "\n",
    "def check_variables(variables, N, sigma, alpha, beta):\n",
    "    try:\n",
    "        y, x_variables, v_variables, w_variables = variables\n",
    "        scaling_variables = pow(10,15)\n",
    "        precision_coefficients = pow(10,20)\n",
    "        if y < 0:\n",
    "            raise InfeasibilityError(\"Negative y variable.\")\n",
    "        for row in v_variables.values():\n",
    "            for v in row.values():\n",
    "                if v < 0:\n",
    "                    raise InfeasibilityError(\"Negative v variable.\")\n",
    "        for row in w_variables.values():\n",
    "            for w in row.values():\n",
    "                if w < 0:\n",
    "                    raise InfeasibilityError(\"Negative w variable.\")\n",
    "        for row in x_variables.values():\n",
    "            for x in row.values():\n",
    "                if x < 0:\n",
    "                    raise InfeasibilityError(\"Negative x variable.\")\n",
    "                    \n",
    "        # checking the dual constraints of type (2)      \n",
    "        y_coefficient_lhs = upper_bound_alpha_plus_delta_1_times_e_to_sigma_beta(N, sigma, alpha, beta)\n",
    "        coefficient_rhs = lower_bound_e_to_minus_one_half_beta(beta) # lower bound on e^(-1/2 beta), scaled by 10^20\n",
    "        for k in range(0, N + 1):\n",
    "            lhs = y_coefficient_lhs * y\n",
    "            for j in range(k, N + 1):\n",
    "                lhs += precision_coefficients * value(v_variables, k, j)\n",
    "            for j in range(0, k + 1):\n",
    "                lhs += precision_coefficients * value(w_variables, j, k)\n",
    "            if k == 1:\n",
    "                lhs += upper_bound_delta_2(N, sigma, alpha, beta) * y # upper bound on delta_2, scaled by 10^20\n",
    "            if k > 0:\n",
    "                for i in range(1, min(k, N -k) + 1):\n",
    "                    lhs += precision_coefficients * value(x_variables, i, k)\n",
    "                for i in range(k, N - k + 1):\n",
    "                    lhs += precision_coefficients * value(x_variables, k, i)\n",
    "                for i in range(1, N + 1):\n",
    "                    j = k - i\n",
    "                    if j >= i and j <= N:\n",
    "                        lhs -= precision_coefficients * value(x_variables, i, j)\n",
    "            rhs = coefficient_rhs * scaling_variables\n",
    "            if lhs > rhs:\n",
    "                print (upper_bound_alpha_plus_delta_1_times_e_to_sigma_beta(N, sigma, alpha, beta))\n",
    "                print (upper_bound_delta_2(N, sigma, alpha, beta))\n",
    "                raise InfeasibilityError(str(lhs) + \" \" + str(rhs) + \" for k=\" + str(k))\n",
    "            y_coefficient_lhs *= upper_bound_e_to_minus_sigma_times_beta(sigma, beta) # upper bound on e^(-sigma * beta), scaled by 10^20\n",
    "            y_coefficient_lhs = y_coefficient_lhs//precision_coefficients + 1\n",
    "            coefficient_rhs *= lower_bound_e_to_minus_beta(beta) # lower bound on e^(- beta), scaled by 10^20\n",
    "            coefficient_rhs = coefficient_rhs//precision_coefficients\n",
    "        \n",
    "        # checking the dual constraints of type (1)   \n",
    "        y_coefficient = upper_bound_4_beta_e_to_sigma_beta(sigma, beta) # upper bound on 4 * beta * e^(sigma * beta), scaled by 10^20\n",
    "        lhs_multiplicator = upper_bound_e_to_minus_sigma_times_beta(sigma, beta) # upper bound on e^(-sigma * beta), scaled by 10^20\n",
    "        for i in range(0, N + 1):\n",
    "            coeff_for_this_i = y_coefficient\n",
    "            for j in range(i, N + 1):\n",
    "                lhs = coeff_for_this_i * y\n",
    "                if (lhs > precision_coefficients * (value(v_variables, i, j) + value(w_variables, i, j))):\n",
    "                    raise InfeasibilityError(\"Inequality for v_ij + w_ij does not hold: lhs \" + str(lhs) + \", rhs \" + str(precision_coefficients *(value(v_variables, i, j) + value(w_variables, i, j))) + \", i: \" + str(i) + \", j: \" + str(j))\n",
    "                coeff_for_this_i *= lhs_multiplicator\n",
    "                coeff_for_this_i = coeff_for_this_i//precision_coefficients + 1\n",
    "            # We need to do the following twice: Once because i increases and once because the next inner loop starts one index later \n",
    "            # Phrased differently, we want y_coefficient to be 4*beta*e^(-(2i-1)*sigma*beta) at the beginning, so from one iteration to the next,\n",
    "            # we need to multiply by e^(-sigma*beta) twice.\n",
    "            y_coefficient *= lhs_multiplicator\n",
    "            y_coefficient = y_coefficient//precision_coefficients + 1\n",
    "            y_coefficient *= lhs_multiplicator\n",
    "            y_coefficient = y_coefficient//precision_coefficients + 1\n",
    "        \n",
    "        objective = lower_bound_sigma_to_minus_2(sigma) * y # the first factor is a lower bound on sigma^(-2), scaled by 10^10\n",
    "        print(\"Solution is feasible and has objective value \" + str(objective) + \"/\" + str(pow(10, 10)*scaling_variables))\n",
    "        upper_bound_sampling = upper_bound_sigma_squared(sigma) * scaling_variables // y + 1 # sigma^2 = 0.439569, so the upper bound is scaled by 10^6*scaling_factor//objective + 1\n",
    "        print(\"This yields an upper bound on the guarantee of the sampling algorithm of less than \" + str(upper_bound_sampling) + \"/10^6\")\n",
    "    except InfeasibilityError as error:\n",
    "        print(\"Solution not feasible; \" + str(error))\n",
    "    except InputError as error:\n",
    "        print(\"InputError; \" + str(error))                "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "486ee304",
   "metadata": {},
   "outputs": [],
   "source": [
    "variables, N, sigma, alpha, beta = read_input(file_path)\n",
    "check_variables(variables, N, sigma, alpha, beta)"
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}