{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "15a8b1c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "using QuantumAlgebra\n",
    "using SymPy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a02f534",
   "metadata": {},
   "source": [
    "### Function definition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "6df38e82",
   "metadata": {},
   "outputs": [],
   "source": [
    "dispnormal(x) = display(\"text/latex\",\"\\$ $(latex(x)) \\\\quad\\\\to\\\\quad $(latex(normal_form(x))) \\$\")\n",
    "acomm(A,B) = A*B + B*A\n",
    "\n",
    "\n",
    "function comm_bo(A::QuantumAlgebra.BaseOperator, B::QuantumAlgebra.BaseOperator)\n",
    "    A_quexpr = QuantumAlgebra.QuExpr(QuantumAlgebra.QuTerm(A))\n",
    "    B_quexpr = QuantumAlgebra.QuExpr(QuantumAlgebra.QuTerm(B))\n",
    "    \n",
    "    c = normal_form(comm(A_quexpr,B_quexpr))   \n",
    "    if c.terms.count==0\n",
    "        c=0\n",
    "    end\n",
    "    return c\n",
    "end\n",
    "\n",
    "\n",
    "function is_selfadjoint(A::QuExpr)\n",
    "    return normal_form(A)==normal_form(adjoint(A))\n",
    "end\n",
    "\n",
    "\n",
    "function CollectTerms(expression::QuExpr)\n",
    "    P = Vector{QuantumAlgebra.Param}[]\n",
    "    B = Vector{QuantumAlgebra.BaseOperator}[]\n",
    "    T = Vector{Int64}()\n",
    "    for (s,t) in expression.terms\n",
    "        if !isempty(s.bares)\n",
    "            push!(B,s.bares.v)\n",
    "            push!(P,s.params)\n",
    "            push!(T,t)\n",
    "        end\n",
    "    end\n",
    "    \n",
    "    BB = [QuantumAlgebra.BaseOpProduct(b) for b in unique(B)]\n",
    "    \n",
    "    Punique = [P[findall(x->x==b.v, B)] for b in BB]\n",
    "    PP = [[QuantumAlgebra.QuTerm(0, QuantumAlgebra.δ[], pint, QuantumAlgebra.ExpVal[], QuantumAlgebra.Corr[],\n",
    "                QuantumAlgebra.BaseOpProduct()) for pint in pext] for pext in Punique]\n",
    "    \n",
    "    TT = [T[findall(x->x==b.v, B)] for b in BB]\n",
    "    return BB, PP, TT\n",
    "end\n",
    "\n",
    "\n",
    "function Terms2QuExpr(BB, PP, TT, separate=false)\n",
    "    AA = Vector{QuantumAlgebra.QuExpr}()\n",
    "    AB = Vector{QuantumAlgebra.QuExpr}()\n",
    "    \n",
    "    for (pext,text) in zip(PP,TT)\n",
    "        A = QuExpr()\n",
    "        for (pint,tint) in zip(pext,text)\n",
    "            QuantumAlgebra._add_with_auto_order!(A,pint,tint)\n",
    "        end\n",
    "        push!(AA,A)\n",
    "    end\n",
    "    \n",
    "    AB = [QuExpr(QuantumAlgebra.QuTerm(0, QuantumAlgebra.δ[], [],\n",
    "                QuantumAlgebra.ExpVal[], QuantumAlgebra.Corr[], b)) for b in BB]\n",
    "    \n",
    "    ### Separation into self-adjoint and npn-selfadjoint parts. Only half of the non-sa part is kept.\n",
    "    if separate \n",
    "        #self-adjoint\n",
    "        i_sa = findall(is_selfadjoint.(AB))\n",
    "        AB_sa = AB[i_sa]\n",
    "        AA_sa = AA[i_sa]\n",
    "        PP_sa = PP[i_sa]\n",
    "        TT_sa = TT[i_sa]\n",
    "        \n",
    "        #not self-adjoint\n",
    "        i_nsa  = deleteat!(collect(range(1,length(AB),step=1)),i_sa)\n",
    "        AB_nsa = AB[i_nsa]\n",
    "        AA_nsa = AA[i_nsa]\n",
    "        PP_nsa = PP[i_nsa]\n",
    "        TT_nsa = TT[i_nsa]\n",
    "        \n",
    "        for i in 1:length(i_nsa)÷2\n",
    "            check = normal_form(adjoint(AB_nsa[i]))\n",
    "            index = findall(x->x==check,AB_nsa)\n",
    "            if check in AB_nsa\n",
    "                deleteat!(AB_nsa, index)\n",
    "            end\n",
    "        end\n",
    "        \n",
    "        keep = findall(in(AB_nsa), AB[i_nsa])\n",
    "        AA_nsa = AA_nsa[keep]\n",
    "        PP_nsa = PP_nsa[keep]\n",
    "        TT_nsa = TT_nsa[keep]        \n",
    "        return AB_sa, AA_sa, PP_sa, TT_sa, AB_nsa, AA_nsa, PP_nsa, TT_nsa \n",
    "    else\n",
    "       return AA, AB \n",
    "    end\n",
    "end\n",
    "\n",
    "\n",
    "function no_exp(Ad, A, n)\n",
    "    F = 0\n",
    "    for i in range(0,n)\n",
    "        F += binomial(n,i) * Ad^(n-i) * A^i\n",
    "    end\n",
    "    return normal_form(F)\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "715b195a",
   "metadata": {},
   "source": [
    "### Parameter initialization \n",
    "\n",
    "All time dependences are included into $\\xi_i$ which should be regarded as\n",
    "$$\n",
    "\\xi_i = e^{i\\omega_i t}\n",
    "$$\n",
    "\n",
    "The goal is to take an expession, expand it, put it in normal form and collect only the non-rotating terms. \n",
    "\n",
    "**TODO** Make the time dependence a method of the bosonic operator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "50198616",
   "metadata": {},
   "outputs": [],
   "source": [
    "ξa = param(:ξ_a,'n')\n",
    "ξb = param(:ξ_b,'n')\n",
    "ξp = param(:ξ_p,'n')\n",
    "\n",
    "@boson_ops b\n",
    "@boson_ops c\n",
    "\n",
    "A  = ξa*a() + ξb*b() + ξp*c()\n",
    "Ad = ξa'*a()'+ ξb'*b()' + ξp'*c()'\n",
    "Q3 = no_exp(Ad, A, 3)\n",
    "Q4 = no_exp(Ad, A, 4)\n",
    "Q5 = no_exp(Ad, A, 5)\n",
    "Q6 = no_exp(Ad, A, 6);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51b833df",
   "metadata": {},
   "source": [
    "### Term collection\n",
    "\n",
    "With `CollectTerms` we put an expression (here $\\phi^2$ or $\\phi^4$) into normal form and collect the different terms. \n",
    "\n",
    "• `BB` contains all different operator compositions.\n",
    "\n",
    "• `PP` for each operator composition in `BB` it takes all its time-dependent prefactors that present in the expression. As the same operator may come with different prefactors due to normal ordering, multiplications, ecc.., `PP` is an array of arrays: one array per each entrance of `BB`.\n",
    "\n",
    "• `TT` numerical prefactor to each entry of `PP`. Has the same structure of `PP`.\n",
    "\n",
    "With `Terms2QuExpr` we separate the self-adjoint components from the non-self-adjoint ones. Because of hermiticity, we can only keep half of the non-self-adjoint components. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "29e01d24",
   "metadata": {},
   "outputs": [],
   "source": [
    "BB, PP, TT = CollectTerms(Q5)\n",
    "AB_sa, AA_sa, PP_sa, TT_sa, AB_nsa, AA_nsa, PP_nsa, TT_nsa  = Terms2QuExpr(BB, PP, TT, true);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89a2872d",
   "metadata": {},
   "source": [
    "### Resonance \n",
    "\n",
    "For each operator composition in `BB` we may have different (one array) time-dependent pre-factors. Below we explore each one individually and check for resonance. \n",
    "\n",
    "We define an array with three entrances, one for each of the frequencies. Each appearance in the prefactor of $\\xi_i$ adds $+1$ in the $i^{th}$, each appearance in the prefactor of $\\xi_i^*$ adds $-1$ in the same.\n",
    "\n",
    "Clearly if all frequencies appear in a \"self-adjoint manner\" their conjugate components cancel the non-conjugate ones leading to `[0,0,0]`: this is the first resonance condition. Example:\n",
    "$$ \\xi_a \\xi_a^* \\,\\,\\xi_b^2\\left(\\xi_b^2\\right)^*  = e^{i(\\omega_a-\\omega_a)t}\\,e^{2i(\\omega_b-\\omega_b)t}\\sim1$$\n",
    "\n",
    "We also have to account for the phase matching condition\n",
    "\n",
    "$$\n",
    "2\\omega_a = \\omega_p + \\omega_b\n",
    "$$\n",
    "\n",
    "This tells us that the `[+2,-1,-1]` and `[-2,+1,+1]` configutations are also resonance configurations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "a2f69123",
   "metadata": {},
   "outputs": [],
   "source": [
    "AB_sa_resonant = Vector{QuantumAlgebra.QuExpr}()\n",
    "AA_sa_resonant = Vector{QuantumAlgebra.QuExpr}()\n",
    "\n",
    "check = fill([], length(PP_sa))\n",
    "for (i,p) in enumerate(PP_sa)\n",
    "    c = []\n",
    "    for pp in p\n",
    "        v = [0,0,0]\n",
    "        for ppp in pp.params\n",
    "            ppp = QuExpr(ppp)\n",
    "            v += [ppp].==[ξa, ξb, ξp]\n",
    "            v -= [ppp].==[ξa', ξb', ξp']\n",
    "        end \n",
    "        append!(c,[v])\n",
    "    end\n",
    "    check[i]=c\n",
    "end\n",
    "\n",
    "\n",
    "for (i,c) in enumerate(check)\n",
    "    A = QuExpr()\n",
    "    success = false\n",
    "    for (j,cc) in enumerate(c)\n",
    "        if cc==[0,0,0] || cc==[2,-1,1] || cc==[-2,1,-1]\n",
    "            success = true\n",
    "            QuantumAlgebra._add_with_auto_order!(A,PP_sa[i][j],TT_sa[i][j])\n",
    "        end\n",
    "    end\n",
    "    if success\n",
    "        push!(AB_sa_resonant, AB_sa[i])\n",
    "        push!(AA_sa_resonant, A)\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "c2655100",
   "metadata": {},
   "outputs": [],
   "source": [
    "AB_nsa_resonant = Vector{QuantumAlgebra.QuExpr}()\n",
    "AA_nsa_resonant = Vector{QuantumAlgebra.QuExpr}()\n",
    "\n",
    "check = fill([], length(PP_nsa))\n",
    "for (i,p) in enumerate(PP_nsa)\n",
    "    c = []\n",
    "    for pp in p\n",
    "        v = [0,0,0]\n",
    "        for ppp in pp.params\n",
    "            ppp = QuExpr(ppp)\n",
    "            v += [ppp].==[ξa, ξb, ξp]\n",
    "            v -= [ppp].==[ξa', ξb', ξp']\n",
    "        end \n",
    "        append!(c,[v])\n",
    "    end\n",
    "    check[i]=c\n",
    "end\n",
    "\n",
    "\n",
    "for (i,c) in enumerate(check)\n",
    "    A = QuExpr()\n",
    "    success = false\n",
    "    for (j,cc) in enumerate(c)\n",
    "        if cc==[0,0,0,0] || cc==[2,-1,1] || cc==[-2,1,-1]\n",
    "            success = true\n",
    "            QuantumAlgebra._add_with_auto_order!(A,PP_nsa[i][j],TT_nsa[i][j])\n",
    "        end\n",
    "    end\n",
    "    if success\n",
    "        push!(AB_nsa_resonant, AB_nsa[i])\n",
    "        push!(AA_nsa_resonant, A)\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "2ccca868",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Resonant self-adjoint terms:\n"
     ]
    }
   ],
   "source": [
    "println(\"Resonant self-adjoint terms:\")\n",
    "\n",
    "for (coeff, oper) in zip(AA_sa_resonant, AB_sa_resonant)\n",
    "    display(\"text/latex\",\"\\$ \\\\Big( $(latex(coeff)) \\\\Big)\\\\,\\\\, $(latex(normal_form(oper))) \\$\")\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "b0a6a4a5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Resonant non-self-adjoint terms (+ h.c):\n"
     ]
    }
   ],
   "source": [
    "println(\"Resonant non-self-adjoint terms (+ h.c):\")\n",
    "\n",
    "for (coeff, oper) in zip(AA_nsa_resonant, AB_nsa_resonant)\n",
    "    display(\"text/latex\",\"\\$ \\\\Big( $(latex(coeff)) \\\\Big)\\\\,\\\\, $(latex(normal_form(oper))) \\$\")\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aeaf2428",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06dd1024",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e216172c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f82c33b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "id": "ab6e7563",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "${a}^\\dagger {b}^\\dagger {a} {b}$"
      ],
      "text/plain": [
       "a†() b†() a() b()"
      ]
     },
     "execution_count": 129,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "collect(AB_sa_resonant[6].terms)[1][1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "id": "b68d3e70",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "${a}^\\dagger {b}^\\dagger {a} {b}$"
      ],
      "text/plain": [
       "a†() b†() a() b()"
      ]
     },
     "execution_count": 130,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v = collect(AB_sa_resonant[6].terms)[1][1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "id": "03df5e3f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4-element Vector{QuantumAlgebra.BaseOperator}:\n",
       " a†()\n",
       " b†()\n",
       " a()\n",
       " b()"
      ]
     },
     "execution_count": 131,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v.bares.v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 217,
   "id": "ff59bc2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "w = v.bares.v\n",
    "\n",
    "if length(w)>1\n",
    "   for i in length(w)÷2+1:length(w)\n",
    "        for j in i-1:-1:1\n",
    "            if comm_bo(w[i],w[j])==0\n",
    "                temp = w[i]\n",
    "                w[i] = w[j]\n",
    "                w[j] = temp\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 218,
   "id": "d59aee79",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4-element Vector{QuantumAlgebra.BaseOperator}:\n",
       " a†()\n",
       " a()\n",
       " b()\n",
       " b†()"
      ]
     },
     "execution_count": 218,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.7.2",
   "language": "julia",
   "name": "julia-1.7"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}