{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example: Standard Quantum Measurement Process Tomography on 1-Qubit with Cvxpy "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from quara.objects.qoperation_typical import (\n",
    "    generate_qoperation,\n",
    "    generate_qoperation_object,\n",
    ")\n",
    "from quara.objects.composite_system_typical import generate_composite_system\n",
    "from quara.protocol.qtomography.standard.standard_qmpt import (\n",
    "    StandardQmpt,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# System\n",
    "num_qubits = 1\n",
    "c_sys = generate_composite_system(mode=\"qubit\", num=num_qubits)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Tester Objects\n",
    "state_names = [\"x0\", \"y0\", \"z0\", \"z1\"]\n",
    "povm_names = [\"x\", \"y\", \"z\"]\n",
    "\n",
    "tester_states = [\n",
    "    generate_qoperation_object(\n",
    "        mode=\"state\", object_name=\"state\", name=name, c_sys=c_sys\n",
    "    )\n",
    "        for name in state_names\n",
    "]\n",
    "tester_povms = [\n",
    "    generate_qoperation_object(\n",
    "        mode=\"povm\", object_name=\"povm\", name=name, c_sys=c_sys\n",
    "    )\n",
    "    for name in povm_names\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# True Object\n",
    "true_object_name = \"z-type2\"\n",
    "true_object = generate_qoperation(\n",
    "    mode=\"mprocess\", name=true_object_name, c_sys=c_sys\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Qmpt\n",
    "on_para_eq_constraint = True\n",
    "qmpt = StandardQmpt(\n",
    "    states=tester_states,\n",
    "    povms=tester_povms,\n",
    "    num_outcomes=true_object.num_outcomes,\n",
    "    on_para_eq_constraint=on_para_eq_constraint,\n",
    "    schedules=\"all\",\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# empi_dists\n",
    "prob_dists = qmpt.calc_prob_dists(true_object)\n",
    "empi_dists = [(10, prob_dist) for prob_dist in prob_dists]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Qmpt with Cvxpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from quara.interface.cvxpy.qtomography.standard.loss_function import (\n",
    "    CvxpyLossFunctionOption,\n",
    "    CvxpyRelativeEntropy,\n",
    "    CvxpyUniformSquaredError,\n",
    "    CvxpyApproximateRelativeEntropyWithZeroProbabilityTerm,\n",
    ")\n",
    "from quara.interface.cvxpy.qtomography.standard.minimization_algorithm import (\n",
    "    CvxpyMinimizationAlgorithm,\n",
    "    CvxpyMinimizationAlgorithmOption,\n",
    ")\n",
    "from quara.interface.cvxpy.qtomography.standard.estimator import (\n",
    "    CvxpyLossMinimizationEstimator,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Loss and Algorithm\n",
    "mode_constraint = \"physical\"\n",
    "#mode_constraint = \"unconstraint\"\n",
    "name_solver = \"mosek\"\n",
    "#name_solver = \"scs\"\n",
    "#name_solver = \"cvxopt\"\n",
    "\n",
    "loss = CvxpyRelativeEntropy()\n",
    "#loss = CvxpyUniformSquaredError()\n",
    "#loss = CvxpyApproximateRelativeEntropyWithZeroProbabilityTerm()\n",
    "loss_option = CvxpyLossFunctionOption()\n",
    "algo = CvxpyMinimizationAlgorithm()\n",
    "algo_option = CvxpyMinimizationAlgorithmOption(\n",
    "    name_solver=name_solver, mode_constraint=mode_constraint\n",
    ")\n",
    "estimator = CvxpyLossMinimizationEstimator()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy.testing as npt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "mode_constraint = physical\n",
      "actual          = [ 5.00000000e-01  4.19112990e-11 -1.71018651e-11  5.00000000e-01\n",
      "  1.82633930e-10 -1.19231459e-11 -4.12192314e-11  1.82634937e-10\n",
      " -9.05795422e-11  4.12194987e-11 -1.19250450e-11 -9.05806765e-11\n",
      "  5.00000000e-01  8.86938019e-11  2.96807284e-11  5.00000000e-01\n",
      " -7.14831436e-12 -7.59382841e-12  2.37046817e-11  7.14824700e-12\n",
      " -5.05432244e-11  2.37044410e-11  7.59362887e-12  5.05433277e-11\n",
      "  5.00000000e-01  4.87125442e-12  6.38844411e-11 -5.00000000e-01]\n",
      "expected        = [ 0.5  0.   0.   0.5  0.   0.   0.   0.   0.   0.   0.   0.   0.5  0.\n",
      "  0.   0.5  0.   0.   0.   0.   0.   0.   0.   0.   0.5  0.   0.  -0.5]\n",
      "squared_error   = 5.490980077466456e-19\n"
     ]
    }
   ],
   "source": [
    "# Estimation\n",
    "sqt = qmpt\n",
    "result = estimator.calc_estimate(\n",
    "    qtomography=sqt,\n",
    "    empi_dists=empi_dists,\n",
    "    loss=loss,\n",
    "    loss_option=loss_option,\n",
    "    algo=algo,\n",
    "    algo_option=algo_option,\n",
    "    is_computation_time_required=True,\n",
    ")\n",
    "var_estimate = result.estimated_var\n",
    "\n",
    "# Test\n",
    "actual = var_estimate\n",
    "expected = true_object.to_var()\n",
    "\n",
    "decimal = 1e-8\n",
    "npt.assert_almost_equal(actual, expected, decimal=decimal)\n",
    "print(\"\\n\")\n",
    "print(\"mode_constraint =\", mode_constraint)\n",
    "print(\"actual          =\", actual)\n",
    "print(\"expected        =\", expected)\n",
    "print(\"squared_error   =\", np.inner(actual-expected, actual-expected))"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "interpreter": {
   "hash": "dcb59d832a18a00f36a6c8652b3e4dd2292d04f0c6e293f19e1cd2c4a4cf8937"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.10"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}