{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example is inspired from a lecture of John Watrous in the [course on Theory of Quantum Information](https://cs.uwaterloo.ca/~watrous/CS766/LectureNotes/08.pdf).\n",
    "\n",
    "The Fidelity between two Hermitian semidefinite matrices P and Q is defined as:\n",
    "\n",
    ">>>> $F(P,Q) = {||{P}^{1/2}{Q}^{1/2} ||}_{tr}$ = max $|trace({P}^{1/2}U{Q}^{1/2})|$\n",
    "\n",
    "where the trace norm $||.||_{tr}$ is the sum of the singular values, and the maximization goes over the set of all unitary matrices U. This quantity can be expressed as the optimal value of the following complex-valued SDP:\n",
    "\n",
    ">>>> **maximize** 1/2 trace(Z+Z*)\n",
    "\n",
    ">>>> subject to $\\left[\\begin{array}{cc}P&Z\\\\{Z}^{*}&Q\\end{array}\\right] \\succeq 0$\n",
    "\n",
    ">>>> where $Z \\in \\mathbf {C}^{n \\times n}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----------------------------------------------------------------------------\n",
      "\tSCS v1.1.8 - Splitting Conic Solver\n",
      "\t(c) Brendan O'Donoghue, Stanford University, 2012-2015\n",
      "----------------------------------------------------------------------------\n",
      "Lin-sys: sparse-direct, nnz in A = 19\n",
      "eps = 1.00e-04, alpha = 1.80, max_iters = 20000, normalize = 1, scale = 5.00\n",
      "Variables n = 9, constraints m = 65\n",
      "Cones:\tprimal zero / dual free vars: 29\n",
      "\tsd vars: 36, sd blks: 1\n",
      "Setup time: 1.45e-04s\n",
      "----------------------------------------------------------------------------\n",
      " Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau | time (s)\n",
      "----------------------------------------------------------------------------\n",
      "     0|      inf       inf      -nan      -inf      -nan       inf  1.39e-04 \n",
      "   100| 2.65e-04  2.64e-02  2.29e-03 -5.19e+00 -5.16e+00  1.37e-15  8.48e-03 \n",
      "   200| 2.07e-06  2.08e-04  1.42e-05 -5.19e+00 -5.19e+00  1.52e-16  1.77e-02 \n",
      "   220| 7.84e-07  7.89e-05  5.38e-06 -5.19e+00 -5.19e+00  2.30e-15  1.96e-02 \n",
      "----------------------------------------------------------------------------\n",
      "Status: Solved\n",
      "Timing: Solve time: 1.96e-02s\n",
      "\tLin-sys: nnz in L factor: 93, avg solve time: 2.19e-06s\n",
      "\tCones: avg projection time: 8.26e-05s\n",
      "----------------------------------------------------------------------------\n",
      "Error metrics:\n",
      "dist(s, K) = 1.1202e-09, dist(y, K*) = 1.7323e-09, s'y/m = -4.9828e-12\n",
      "|Ax + s - b|_2 / (1 + |b|_2) = 7.8427e-07\n",
      "|A'y + c|_2 / (1 + |c|_2) = 7.8923e-05\n",
      "|c'x + b'y| / (1 + |c'x| + |b'y|) = 5.3752e-06\n",
      "----------------------------------------------------------------------------\n",
      "c'x = -5.1910, -b'y = -5.1910\n",
      "============================================================================\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "5.191031665238333"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Convex\n",
    "n = 20\n",
    "P = randn(n,n) + im*randn(n,n)\n",
    "P = P*P'\n",
    "Q = randn(n,n) + im*randn(n,n)\n",
    "Q = Q*Q'\n",
    "Z = ComplexVariable(n,n)\n",
    "objective = 0.5*real(trace(Z+Z'))\n",
    "constraint = [P Z;Z' Q] ⪰ 0\n",
    "problem = maximize(objective,constraint)\n",
    "solve!(problem)\n",
    "computed_fidelity = evaluate(objective)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×2 Array{Complex{Float64},2}:\n",
       " 1.73143-5.55112e-17im        -0.047971+0.182073im\n",
       "        -0.047971-0.182073im   0.951632+0.0im     "
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Verify that computer fidelity is equal to actual fidelity\n",
    "P1,P2 = eig(P)\n",
    "sqP = P2 * diagm([p1^0.5 for p1 in P1]) * P2'\n",
    "Q1,Q2 = eig(Q)\n",
    "sqQ = Q2 * diagm([q1^0.5 for q1 in Q1]) * Q2'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.1910314597465606"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "actual_fidelity = sum(svd(sqP * sqQ)[2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the actual fidelity value is very close the computed fidelity value."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Julia 0.5.0-rc3",
   "language": "julia",
   "name": "julia-0.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}