Non_Diagonal_Shifts.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import division, print_function\n",
    "from sympy import Rational, re, im\n",
    "import math, cmath, time\n",
    "\n",
    "from CFT import *\n",
    "from Correlators import *\n",
    "from Spectrum_bootstrap import *\n",
    "from Data import Table\n",
    "import numbers\n",
    "from Auxiliary_classes import numerical\n",
    "import numbers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Shift equations and four-point functions in non-diagonal CFT\n",
    "In this notebook we define classes to compute four-point functions and structure constants in a non-diagonal CFT.\n",
    "\n",
    "The structure constants are computed by recursively using the shift equations of the analytic conformal bootstrap.\n",
    "\n",
    "Four-point functions are computed using the non-diagonal spectrum $\\mathbb{S}_{2\\mathbb{Z}, \\mathbb{Z}+\\frac{1}{2}}$. \n",
    "\n",
    "##  Notations\n",
    "\n",
    "We write the central charge in terms of a parameter $\\beta$ such that\n",
    "\\begin{align}\n",
    "c = 1-6\\left(\\beta - \\frac{1}{\\beta}\\right)^2\\, .\n",
    "\\end{align}\n",
    "Diagonal fields $V^D_{P}$ can have arbitrary mommentums $P=\\bar P$.\n",
    "Non-diagonal fields $V^N_{(r,s)}$ are labelled by two half-integer numbers $r,s$, and their left and right momentums are\n",
    "\\begin{align}\n",
    "P &= P_{(r,s)} = \\frac{1}{2}\\left(r\\beta - \\frac{s}{\\beta}\\right)\\, ,\\\\\n",
    "\\bar{P} &= P_{(r,-s)}\\, .\n",
    "\\end{align}\n",
    "\n",
    "The central charge is assumed to obey $\\Re{c}<13$.\n",
    "\n",
    "##  Shift equations and four-point structure constants\n",
    "\n",
    "The ratios that appear in shift equations corresponding to the degenerate field are $V_{\\langle 2,1 \\rangle}$ are:\n",
    "\\begin{align}\n",
    "\\rho\\left(V_1|V_2, V_3 \\right) = -(-1)^{2s_2} \\frac{\\Gamma(-2\\beta P_1)\\Gamma(-2\\beta \\sigma_1\\bar{P}_1)}{\\Gamma(2\\beta P_1)\\Gamma(2\\beta \\sigma_1\\bar{P}_1)} \\prod_{\\pm \\pm}\\frac{\\Gamma(\\frac{1}{2} + \\beta P_1 \\pm \\beta P_2 \\pm \\beta P_3)}{(\\frac{1}{2} - \\beta \\sigma_1 \\bar{P}_1 \\pm \\beta \\bar{P}_2 \\pm \\beta \\bar{P}_3)}\\, ,\n",
    "\\end{align}\n",
    "and \n",
    "\\begin{align}\n",
    "\\rho(V) = -\\frac{\\Gamma(-2\\beta P)\\Gamma(-2\\beta \\sigma\\bar{P})}{\\Gamma(2\\beta P)\\Gamma(2\\beta \\sigma\\bar{P})}  \\frac{\\Gamma(\\beta^2+2\\beta P)\\Gamma(1-\\beta^2+2\\beta P)}{\\Gamma(\\beta^2-2\\beta \\sigma \\bar{P})\\Gamma(1-\\beta^2-2\\beta \\sigma \\bar{P})}\\, .\n",
    "\\end{align}\n",
    "\n",
    "The shift equations for the four-point structure constants are\n",
    "\\begin{align}\n",
    "\\frac{D(V^{+}_s)}{D(V^{-}_s)} = \\frac{\\rho \\left(V_s|V_1,V_2 \\right)\\rho \\left(V_s|V_3,V_4 \\right)}{\\rho\\left( V_s \\right)}\\, .\n",
    "\\end{align}\n",
    "\n",
    "We also implement the dual shift equations, corresponding to the degenerate field $V_{\\langle 1,2 \\rangle}$.\n",
    "\n",
    "\n",
    "## Normalizations\n",
    "\n",
    "Computing structure constants using the shift equations makes it necessary to choose a starting point for the recursions. We can arbitrarily fix structure constants corresponding to the lowest dimension to be $1$, i.e:\n",
    "\n",
    "\\begin{align}\n",
    "D\\left( V^N_{\\left( 0,\\frac{1}{2} \\right)} \\right) = 1\\, .\n",
    "\\end{align}\n",
    "\n",
    "However, comparing $s$-and $t$-channel computations requires that we determine the relative normalization between both channels. This is given by\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{D_{2341}\\left( V^N_{\\left( 0,\\frac{1}{2} \\right)} \\right)}{D_{1234}\\left( V^N_{\\left( 0,\\frac{1}{2} \\right)} \\right)} = (-1)^{\\sum_{i = 1}^{4}\\text{Spin}_i } \\prod_{j=0}^{m-1} \\prod_{k=0}^{n-1} \\left[ \\frac{\\tilde{\\rho} \\left( V^N_{\\left( r_2, \\, s_2 -\\delta_s (2j+1) \\right) }, V^N_{\\left( 0,\\frac{1}{2} \\right)} , V_3 \\right) } {\\tilde{\\rho}\\left( V^N_{\\left( r_2, \\, s_2 -\\delta_s (2j+1) \\right) }, V^N_{\\left( 0,\\frac{1}{2} \\right)}, V_1 \\right) }\\right]^{\\delta_s}\n",
    "\\times\n",
    "\\left[ \\frac{\\rho \\left( V^N_{\\left( r_2 -\\delta_r (2k+1) , \\, s_4 \\right) }, V^N_{\\left( 0,\\frac{1}{2} \\right)}, V_3 \\right) } {\\rho\\left( V^N_{\\left( r_2 -\\delta_r (2k+1) , \\, s_4 \\right) }, V^N_{\\left( 0,\\frac{1}{2} \\right)}, V_1 \\right) }\\right]^{\\delta_r} \\, , \n",
    "\\end{align}\n",
    "\n",
    "where \n",
    "\n",
    "\\begin{align}\n",
    "s_2 - s_4 &= 2 \\delta_s m, \\quad m\\in \\mathbb{Z}_{\\geq 0}, \\\\\n",
    "r_2 - r_4 &= 2 \\delta_r n, \\quad n\\in \\mathbb{Z}_{\\geq 0}\\, .\n",
    "\\end{align}\n",
    "\n",
    "Thus, when computing four-point functions we normalize the first $s$-channel structure constant as $D_{1234}\\left( V^N_{\\left( 0,\\frac{1}{2} \\right)} \\right) = 1$, and multiply the $t$-channel computation by the factor above.\n",
    "\n",
    "### Ordering of the fields\n",
    "\n",
    "Four-point functions will be computed correctly only if diagonal and non-diagonal fieds alternate, i.e. if they are of the types\n",
    "\\begin{align}\n",
    "\\left\\langle  V^D_1 V^N_2 V^D_3 V^N_4 \\right\\rangle \\, ,\\quad \\left\\langle  V^N_1 V^D_2 V^N_3 V^D_4 \\right\\rangle\\, .\n",
    "\\end{align}\n",
    "\n",
    "These are the cases where the $s$- and $t$-channel spectrums are both $\\mathbb{S}_{2\\mathbb{Z}, \\mathbb{Z}+\\frac{1}{2}}$.\n",
    "\n",
    "## Four-point structure constants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class ND_csts:\n",
    "    \"\"\" Computing the structure constants of a given four point function.  \"\"\"\n",
    "    \n",
    "    def __init__(self, ext_fields):\n",
    "        \"\"\"ext_fields: LIST of four FIELDS. All should have the same central charge.\n",
    "           The value of the central charge is taken form the first field. \"\"\" \n",
    "        \n",
    "        self.fields = ext_fields\n",
    "        \n",
    "        self.Ps = [field.dims[0].get('P') for field in ext_fields ]\n",
    "        self.bPs = [field.dims[1].get('P') for field in ext_fields]\n",
    "                      \n",
    "        self.charge = self.fields[0].dims[0].charge\n",
    "        self.beta = self.charge.get('beta')\n",
    "        \n",
    "    @staticmethod\n",
    "    def rho(Vs, dual = False, full = True):\n",
    "        \"\"\"Compute rho and rho-tilde.\n",
    "        Vs: LIST containing three or a single field (for rho(V)).\n",
    "        dual = If True, changes beta --> -1/beta in the arguments of the gammas, giving rho-tilde.\n",
    "        deg = If True, generates rho(V). \n",
    "        full = If True, computes prefectors de\n",
    "        Use: By default, give three fields and get rho(V1|V2,V3).\n",
    "        dual = True --> get rho-tilde(V1|V2,V3).\n",
    "        Degenerate rho's: Set deg = True, and give only one field. The\n",
    "        degenerate momemtums are created within rho. \"\"\"\n",
    "        \n",
    "        \n",
    "        beta = Vs[0].dims[0].charge.get('beta') # We don't use self.beta because it's a static method.\n",
    "        \n",
    "        deg = (len(Vs) == 1)\n",
    "               \n",
    "        sigma = 1\n",
    "        \n",
    "        if dual:\n",
    "            \"\"\"Inverts beta to get rho tilde. If the shifted field is non-diagonal,\n",
    "            changes the sign of sigma.\n",
    "            We determine if the field is diagonal by comparing the left and right Dimension objects\n",
    "            (not the value, the instance of the class.)\"\"\"\n",
    "            beta = -1/beta\n",
    "            if Vs[0].non_diagonal:\n",
    "                sigma = -1\n",
    "        \n",
    "        if deg:\n",
    "            Pdeg = Dimension(variable = 'degenerate', value = (2,1), charge = Charge('beta', beta )).get('P')\n",
    "            \n",
    "            L_P = [Vs[0].dims[0].get('P'), Pdeg, Vs[0].dims[0].get('P')]\n",
    "            R_P = [Vs[0].dims[1].get('P'), Pdeg, Vs[0].dims[1].get('P')]\n",
    "        else:\n",
    "            L_P = [V.dims[0].get('P') for V in Vs]\n",
    "            R_P = [V.dims[1].get('P') for V in Vs]\n",
    "        \n",
    "        \n",
    "        \"\"\"Performing the computation. This takes a few steps:\n",
    "        1 - Compute the factor depending on all momentums\n",
    "        2 - if full\n",
    "        Determine the sign factor.\n",
    "        Compute the factor depending only on the shifted field.\n",
    "        3 - return product of computed factors.\n",
    "        \"\"\"\n",
    "             \n",
    "        prod = 1\n",
    "        for i in [+1,-1]:\n",
    "            for j in [+1,-1]:\n",
    "                prod *= (scipy.special.gamma(1/2 + beta* L_P[0] + i * beta*L_P[1]+ j * beta*L_P[2])\n",
    "                         /scipy.special.gamma(1/2 - sigma *beta* R_P[0] + i * beta*R_P[1]+ j * beta * R_P[2]))\n",
    "                  \n",
    "        if full:\n",
    "            \"\"\"Computes sign and normalization factor\"\"\"\n",
    "            if deg:\n",
    "                sign =1\n",
    "            elif dual:\n",
    "                sign = Vs[1].signs[0]\n",
    "            else:\n",
    "                sign = Vs[1].signs[1]\n",
    "                \n",
    "            norm = - (scipy.special.gamma(-2*beta*L_P[0])*scipy.special.gamma(-2*beta*sigma *R_P[0]) \n",
    "                   /scipy.special.gamma(2*beta*L_P[0]) /scipy.special.gamma(2*beta*sigma *R_P[0]) )\n",
    "        \n",
    "            return sign * norm * prod\n",
    "        else:\n",
    "            return prod\n",
    "\n",
    "    def relative_norm(self):\n",
    "        \"\"\"Relative normalization between s and t channels.\"\"\"\n",
    "        \n",
    "        D0 = Dimension(variable = 'degenerate', value = (0,1/2), charge = Charge('beta', self.beta))\n",
    "        V0 = Field(dim = D0, non_diagonal = True)\n",
    "        \n",
    "        nd_fields = [] # a list of the non-diag fields\n",
    "        d_fields = [] # a list of the diag fields\n",
    "        for f in self.fields:\n",
    "            if f.non_diagonal == True:\n",
    "                nd_fields.append(f)\n",
    "            else:\n",
    "                d_fields.append(f)\n",
    "        \n",
    "                        \n",
    "        (r2, s2) = nd_fields[0].dims[0].get('degenerate') \n",
    "        (r4, s4) = nd_fields[-1].dims[0].get('degenerate')\n",
    "        \n",
    "                \n",
    "        \"\"\"Testing if a reflection is necessary for the shifts in s.\"\"\"\n",
    "        if (s2-s4) % 2 != 0:\n",
    "            (r4, s4) = (-r4, -s4)\n",
    "                \n",
    "        dr = 1 if r2-r4 >= 0 else -1\n",
    "        ds = 1 if s2-s4 >= 0 else -1\n",
    "        \n",
    "        m = int(abs(s2-s4)/2)\n",
    "        n = int(abs(r2-r4)/2)\n",
    "        \n",
    "        s_factor = 1\n",
    "        r_factor = 1\n",
    "        \n",
    "        for j in range(m):\n",
    "            sj = s2 - ds*(2*j+1)\n",
    "            \n",
    "            Dj = Dimension(variable = 'degenerate', value = (r2, sj), charge = Charge('beta', self.beta))\n",
    "            Vj = Field(dim = Dj, non_diagonal = True)\n",
    "            \n",
    "            s_factor *= (ND_csts.rho([Vj, V0, d_fields[-1]], dual = True, full = False)\n",
    "                         /ND_csts.rho([Vj, V0, d_fields[0]], dual = True, full = False) )**ds\n",
    "        \n",
    "        for k in range(n):\n",
    "            rk = r2 - dr*(2*k+1)\n",
    "            \n",
    "            Dk = Dimension(variable = 'degenerate', value = (rk, s4), charge = Charge('beta', self.beta))\n",
    "            Vk = Field(dim = Dk, non_diagonal = True)\n",
    "            \n",
    "            r_factor *= (ND_csts.rho([Vk, V0, d_fields[-1]], full = False)\n",
    "                         /ND_csts.rho([Vk, V0, d_fields[0]], full = False) )**dr\n",
    "        \n",
    "        spin = nd_fields[0].spin + nd_fields[-1].spin\n",
    "        \n",
    "        return numerical( (-1)**spin *r_factor * s_factor )\n",
    "    \n",
    "    def Shift(self, pair, variable = 'r'):\n",
    "        \"\"\"Takes the indices of a field, and a variable as a string. \n",
    "        Returns D(r+1,s)/D(r-1,s) or D(r,s+1)/D(r,s-1).\"\"\"\n",
    "        \n",
    "        Drs = Dimension(variable = 'degenerate', value = pair, charge = Charge('beta', self.beta))\n",
    "        Vrs = Field(dim = Drs, non_diagonal = True)\n",
    "        \n",
    "        if variable == 'r':                         \n",
    "            return (ND_csts.rho([Vrs, self.fields[0], self.fields[1]])\n",
    "                    *ND_csts.rho([Vrs, self.fields[2], self.fields[3]])\n",
    "                   /ND_csts.rho([Vrs]) )\n",
    "        \n",
    "        if variable == 's':    \n",
    "            return (ND_csts.rho([Vrs, self.fields[0], self.fields[1]], dual = True)\n",
    "                    *ND_csts.rho([Vrs, self.fields[2], self.fields[3]], dual = True)\n",
    "                   /ND_csts.rho([Vrs], dual = True) )\n",
    "        \n",
    "    def value(self, pair, path = True):\n",
    "        \"\"\"Calculating the D(r,s)/D(0,1/2), for a single choice of indices.\n",
    "        If path is True, first the shift in r is done. If False, first we do s.\"\"\"\n",
    "        \n",
    "        cst = 1\n",
    "        (r, s) = pair\n",
    "        \n",
    "        sign_r = 1 if r >= 0 else -1\n",
    "        sign_s = 1 if s >= 0 else -1\n",
    "        \n",
    "        s0 = 1/2 if (s-1/2) % 2 == 0 else -1/2\n",
    "        s_start = int(s0 + sign_s/2)\n",
    "        s_end = int(s -sign_s/2)\n",
    "                       \n",
    "        if path == True:\n",
    "            for rstep in range(sign_r*1, r,sign_r*2):\n",
    "                cst *=self.Shift((rstep, s0), 'r')**sign_r\n",
    "\n",
    "            for sstep in range(s_start, s_end, sign_s*2):\n",
    "                cst *=self.Shift((r, sstep + sign_s/2), 's')**sign_s\n",
    "        else:\n",
    "            for sstep in range(s_start, s_end, sign_s*2):\n",
    "                cst *=self.Shift((0, sstep + sign_s/2), 's')**sign_s\n",
    "            for rstep in range(sign_r*1, r ,sign_r*2):\n",
    "                cst *=self.Shift((rstep, s), 'r')**sign_r\n",
    "            \n",
    "        return cst\n",
    "    \n",
    "    def make_table(self, bounds, tex = False, tex_cmd = False, real = True):\n",
    "        \"\"\" Takes (positive) bounds for (r,s) and builds the table of structure\n",
    "        constants. It also gives two lists with all the values of r and s, to use\n",
    "        as reference when looking for constants in the table.\n",
    "        tex = If True, builds a tex table and displays it.\n",
    "        tex_cmd = Prints the Tex command.\n",
    "        real = Keeps only the real values of the structure constants.\"\"\"\n",
    "        \n",
    "        (rbound,sbound) = bounds\n",
    "        s_end = int(sbound +1/2)\n",
    "        \n",
    "        table = []\n",
    "        rowp = [1] # Starting point (0,1/2)\n",
    "        rowm = [1] # Starting point (0,-1/2)\n",
    "        \n",
    "        r_ref = [0] #Initial value for r indices.\n",
    "        s_ref = [-1/2,1/2] # Initial values for s indices.\n",
    "        for r in range(1, rbound, 2):\n",
    "            \n",
    "            rowm.insert(0, rowm[0]/self.Shift((-r,-1/2), 'r'))\n",
    "            rowm.append(rowm[-1]*self.Shift((r,-1/2), 'r'))\n",
    "                        \n",
    "            rowp.insert(0, rowp[0]/self.Shift((-r,1/2), 'r'))\n",
    "            rowp.append(rowp[-1]*self.Shift((r,1/2), 'r'))\n",
    "            \n",
    "            r_ref.insert(0,-r-1)\n",
    "            r_ref.append(r+1)\n",
    "            \n",
    "        table.extend([rowm, rowp])\n",
    "                \n",
    "        for s in range(1, s_end):\n",
    "            \n",
    "            new_row_up = []\n",
    "            new_row_down = []\n",
    "            \n",
    "            for j in range(0,rbound+1):\n",
    "                new_row_down.append(table[1][j]/self.Shift((-rbound+2*j,-s+1/2 ), 's'))\n",
    "                new_row_up.append(table[-2][j]*self.Shift((-rbound+2*j,s-1/2), 's') )\n",
    "                \n",
    "                            \n",
    "            table.insert(0,new_row_down)\n",
    "            table.append(new_row_up)\n",
    "            \n",
    "            s_ref.insert(0, -s-1/2)\n",
    "            s_ref.append(s+1/2)\n",
    "        \n",
    "        \"\"\"making the tex table \"\"\"\n",
    "        if tex or tex_cmd:\n",
    "            \n",
    "            pretty_csts = []\n",
    "            for row in table:\n",
    "                pretty_row = []\n",
    "                for element in row:\n",
    "                    if abs(element) < 10**(-4):\n",
    "                        pretty_row.append(Table.format_exp('{:.8g}'.format(element)))\n",
    "                    else:\n",
    "                        pretty_row.append(Table.format_result(element, digits = 12, real = real) )\n",
    "                pretty_csts.append(pretty_row)\n",
    "            \n",
    "            tbl_st = r'\\begin{array}{|c|'\n",
    "            \n",
    "            for i  in r_ref:\n",
    "                tbl_st += r'c'\n",
    "\n",
    "            tbl_st += r'| } \\hline '+ '\\n' + r's \\backslash r ' + '\\n'\n",
    "            \n",
    "            for r in r_ref:\n",
    "                tbl_st += ' & {: d}'.format(r)\n",
    "            \n",
    "            tbl_st += r'\\\\ \\hline ' + '\\n'\n",
    "            \n",
    "            for i in range(len(s_ref)):\n",
    "                tbl_st += latex(Rational(s_ref[i])) \n",
    "                for j in range(len(r_ref)):\n",
    "                    tbl_st += r' &' + pretty_csts[i][j]\n",
    "                \n",
    "                tbl_st += r'\\\\ ' + '\\n'\n",
    "            tbl_st += r'\\hline ' + '\\n' + r'\\end{array}'\n",
    "            \n",
    "            if tex:\n",
    "                display(Math(tbl_st))\n",
    "            \n",
    "            if tex_cmd:\n",
    "                print( tbl_st )\n",
    "            \n",
    "        return [table,r_ref,s_ref]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Four-point functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class ND4:\n",
    "    \"\"\"Calculates 4 point functions with the spectrum (2Z, Z+1/2) in the s- or t-channel.\"\"\"\n",
    "    \n",
    "    def __init__(self, ext_fields, L = 5, blocklevel = 12):\n",
    "        \"\"\"ext_fields: LIST of four Fields. \n",
    "        L = Truncation parmeter of the spectrum\n",
    "        blocklevel = Truncation paramenter of the conformal blocks.\n",
    "        \"\"\"\n",
    "        \n",
    "        self.beta = ext_fields[0].dims[0].charge.get('beta')\n",
    "        self.ext_fields = ext_fields\n",
    "        self.L = L\n",
    "        self.blocklevel = blocklevel\n",
    "        \n",
    "        \"\"\"Building the internal spectrum and it's truncation. \"\"\"    \n",
    "        self.spectrum = Another_Spectrum(self.beta)\n",
    "        self.int_fields = self.spectrum.truncate(self.L)\n",
    "        \n",
    "        \"\"\"Building blocks\"\"\"\n",
    "        dims_l = [f.dims[0] for f in self.ext_fields] # Left dims\n",
    "        dims_r = [f.dims[1] for f in self.ext_fields] # Right dims\n",
    "        \n",
    "        \"\"\"s-blocks\"\"\"\n",
    "        self.s_blocks = [Block(Nmax = self.blocklevel, dims = dims_l),\n",
    "                        Block(Nmax = self.blocklevel, dims = dims_r)]\n",
    "        \"\"\"t-blocks\"\"\"\n",
    "        self.t_blocks = [Block(t_channel = True, Nmax = self.blocklevel, dims = dims_l),\n",
    "                        Block(t_channel = True, Nmax = self.blocklevel, dims = dims_r)]\n",
    "                \n",
    "        \"\"\"Finding maximal indices fot the constants table\"\"\"\n",
    "            \n",
    "        r_vals = [f.dims[0].get('degenerate')[0] for f in self.int_fields]\n",
    "        s_vals = [f.dims[0].get('degenerate')[1] for f in self.int_fields]\n",
    "        \n",
    "        self.rmax = max(r_vals)\n",
    "        self.smax = max(s_vals)\n",
    "            \n",
    "        self.csts = self.constants() #Table of s- and t-ch constants, and normalization.\n",
    "    \n",
    "    def constants(self):\n",
    "        \"\"\"Makes a table with the constants, and references.\n",
    "        Computes the relative normalization between s and t channel\"\"\"\n",
    "        \n",
    "        \"\"\"s-channel constants and normalization\"\"\"\n",
    "        s_csts = ND_csts(self.ext_fields)\n",
    "        \n",
    "        N = s_csts.relative_norm() # Relative normalization\n",
    "        s_table = s_csts.make_table((self.rmax, self.smax)) #table of s-constants\n",
    "        \n",
    "        \"\"\"t-channel constants\"\"\"\n",
    "        \n",
    "        t_csts = ND_csts([self.ext_fields[1],self.ext_fields[2],\n",
    "                        self.ext_fields[3], self.ext_fields[0]])\n",
    "        # The permutation chosen for t corresponds to the one on the Block class. \n",
    "        \n",
    "        t_table = t_csts.make_table((self.rmax, self.smax)) #table of t-constants\n",
    "        \n",
    "        \n",
    "        return [s_table[0],t_table[0], s_table[1],s_table[2], N] # s-constants, t-constants, ref_r, ref_s\n",
    "  \n",
    "    def blocks(self, z, reduced = True, channel = None):\n",
    "        \"\"\" Computes the needed values of conformal blocks. Computes both channels\n",
    "        if None, only the s-channel if True, and only the t-channel if False.\n",
    "        \"\"\"\n",
    "        \n",
    "        q = mpmath.qfrom(m = z)\n",
    "        if channel or channel is None:\n",
    "            block_s_l = BlockNum(block = self.s_blocks[0], q = q)\n",
    "            block_s_r = BlockNum(block = self.s_blocks[1], q = q.conjugate())\n",
    "        if not channel or channel is None:\n",
    "            block_t_l = BlockNum(block = self.t_blocks[0], q = q)\n",
    "            block_t_r = BlockNum(block = self.t_blocks[1], q = q.conjugate())\n",
    "        values = []\n",
    "      \n",
    "        for field in self.int_fields:\n",
    "            value = [field.dims[0].get('degenerate')] #Reference indices.\n",
    "            if channel or channel is None:\n",
    "                value.append(block_s_l.value(field.dims[0], reduced) \n",
    "                         * block_s_r.value(field.dims[1], reduced))\n",
    "            if not channel or channel is None:\n",
    "                value.append(block_t_l.value(field.dims[0], reduced)\n",
    "                             *block_t_r.value(field.dims[1], reduced))\n",
    "            \n",
    "            values.append(value)\n",
    "            \n",
    "        return values\n",
    "\n",
    "    def values(self, z, channel = None, norm = True):\n",
    "                \n",
    "        blocks = self.blocks(z = z, channel = channel) #List of indices and block values.\n",
    "        \n",
    "        table = self.csts[:4] # Tables of constants and references\n",
    "        \n",
    "        N = self.csts[4] # t-channel normalization\n",
    "        \n",
    "        channels = [0, 1] if channel is None else [0]\n",
    "        values = [0 for i in channels]\n",
    "                               \n",
    "        for data in blocks:\n",
    "            \"\"\"Determine the position of the constant in the table. Same for both channels.\"\"\"\n",
    "            index_r = table[2].index(data[0][0]) #Determining r index\n",
    "            index_s = table[3].index(data[0][1]) #Determining s index\n",
    "            \n",
    "            if channel or channel is None:\n",
    "                values[0] += table[0][index_s][index_r] * data[1] #s-channel values\n",
    "            if not channel or channel is None:\n",
    "                values[-1] += table[1][index_s][index_r] * data[-1] #t-channel values\n",
    "        \n",
    "        if norm:\n",
    "            if not channel or channel is None:\n",
    "                values[-1] *= N # multiplying by the t-s normalization\n",
    "        \n",
    "        if channel is not None:\n",
    "            values = values[0] #If only one channel, returns a number instead of a list.\n",
    "                        \n",
    "        return values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  An auxiliary class for building arguments of four-point functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class four_point_args:\n",
    "    \"\"\"A class for building arguments of four-point functions. \n",
    "    Takes a central charge and a list of fields or paramenters to build them, and creates a \n",
    "    list of fields. The method channel allows to permute them.\n",
    "    \"\"\"\n",
    "          \n",
    "    def __init__(self, beta, *args):\n",
    "        \"\"\" The first argument is the parameter Beta of the central charge. The rest can be:\n",
    "        Field: Just add it to a list.\n",
    "        Number: Creates a diagonal field with this value of Delta, and the given central charge.\n",
    "        Adds it to the list.\n",
    "        Tuple (r,s, string): (r,s) are the indices of a field. The string should be 'd' or 'n'\n",
    "        specifying if the field created is diagonal or not. Adds the field to the list.\n",
    "        \n",
    "        If given two arguments, it duplicates them to create a list of two fields of the type\n",
    "        1212.\n",
    "        The created list has the fields ordered in the same way as the arguments. \n",
    "        \"\"\"\n",
    "        self.args = args\n",
    "        self.beta = beta\n",
    "        self.charge = Charge('beta', beta)\n",
    "        \n",
    "        if len(args) not in [2,4]:\n",
    "            print( 'WARNING: Wrong number of arguments' )\n",
    "                    \n",
    "        fields = []\n",
    "        \n",
    "        for a in self.args:\n",
    "            \n",
    "            if isinstance(a, Field):\n",
    "                \"\"\"If it gets a Field, it adds it to the list.\"\"\"\n",
    "                fields.append(a)\n",
    "            elif isinstance(a, numbers.Number):\n",
    "                \"\"\"If it gets a Number, creates a diagonal field of this Delta.\"\"\"\n",
    "                d = Dimension(variable = 'Delta', value = a, charge = self.charge)\n",
    "                fields.append(Field(dim = d, diagonal = True))\n",
    "            elif isinstance(a, tuple):\n",
    "                \"\"\"If it gets a tuple, it creates a diag or non-diag field with the first\n",
    "                 two values of the tuple as indices\"\"\"\n",
    "                d = Dimension(variable = 'degenerate', value = (a[0],a[1]), charge = self.charge)\n",
    "                if a[2] =='d':\n",
    "                    fields.append(Field(dim = d, diagonal = True))\n",
    "                elif a[2] =='n':\n",
    "                    fields.append(Field(dim = d, non_diagonal = True))\n",
    "        \n",
    "        if len(fields) == 2:\n",
    "            fields.extend(fields)\n",
    "        \n",
    "        self.fields = fields\n",
    "        \n",
    "    @staticmethod\n",
    "    def permute(sigma, source):\n",
    "        \"\"\"Takes a string of numbers between 1 and 4 and a list, and gives a new list\n",
    "        with elements permuted.\n",
    "        sigma : A STRING of the type 1234 or 4132 or any permutation. \n",
    "        source: The LIST whose elements we want to permute \"\"\"\n",
    "\n",
    "        permuted = []\n",
    "        for i in list(sigma):\n",
    "            permuted.append(source[int(i)-1])\n",
    "        return permuted\n",
    "      \n",
    "    def channel(self, perm):\n",
    "        \"\"\"Takes a string as an argument, and returns the list of \n",
    "        fields according to it. The string can be 's', 't', 'u' or\n",
    "        any permutation of '1234'.\n",
    "        \"\"\"\n",
    "        channels = {'s' : '1234', 't': '2341', 'u': '1342'}\n",
    "        # The permutation chosen for t corresponds to the one on the Block class. \n",
    "        \n",
    "        if perm in channels:\n",
    "            return four_point_args.permute(channels[perm], self.fields)\n",
    "        else:\n",
    "            return four_point_args.permute(perm, self.fields)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tests\n",
    "\n",
    "\n",
    "## Recovering numerical bootstrap results\n",
    "\n",
    "We compute four-point structure constants in the case of $c=0$ and four $(0,\\frac12)$ fields (two diagonal, two non-diagonal). These can be compared with earlier numerical bootstrap results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if __name__ == '__main__':\n",
    "    \n",
    "    \"\"\"Computing the values of certain structure constants. Results using the methods 'value'\n",
    "    and 'make_table' of the class ND_csts.\n",
    "    \"\"\"\n",
    "    c = 0\n",
    "    C = Charge('c', c)\n",
    "    beta0 = C.get('beta')\n",
    "    \n",
    "    Champs = four_point_args(beta0, (0,1/2,'d'), (0,1/2,'n'))\n",
    "    \n",
    "    pairs = [(0,1/2), (2,1/2), (0,3/2), (2,3/2), (0,5/2), (4,1/2), (2,5/2), (4,3/2), (0,7/2)]\n",
    "    \n",
    "    Csts = ND_csts(Champs.channel('s'))\n",
    "    \n",
    "    \"\"\"Results using the method 'value' . \"\"\"\n",
    "    for pair in pairs:\n",
    "        print( str(pair), Csts.value(pair) )\n",
    "\n",
    "    \n",
    "    \"\"\"Table in tex format\"\"\"\n",
    "    Csts.make_table((4,5/2), tex = True, tex_cmd = True ) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Crossing symmetry of four-point functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if __name__ == '__main__':\n",
    "    \n",
    "    \"\"\"Case of c=0\"\"\"\n",
    "    \n",
    "    c = 0\n",
    "    eps = 10**(-9)\n",
    "    C = Charge('c', c+eps)\n",
    "    beta0 = C.get('beta')\n",
    "    \n",
    "    Champs = four_point_args(beta0, (0,1/2,'d'), (0,1/2,'n')).channel('s')\n",
    "    \n",
    "    F = ND4(Champs, L = 8, blocklevel = 12)\n",
    "    \n",
    "    for x in range(1,5):\n",
    "        val_x = F.values(x/10, channel = None)\n",
    "        \n",
    "        print( 'x = %s' % (x/10) )\n",
    "        print( 's-ch: %s' % val_x[0] )\n",
    "        print( 't-ch: %s' % val_x[1] )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if __name__ == '__main__':\n",
    "    \n",
    "    \"\"\"Generic charge and dimensions\"\"\"\n",
    "    \n",
    "    c = 4.72+0.12*1j\n",
    "    eps = 0\n",
    "    C = Charge('c', c+eps)\n",
    "    beta0 = C.get('beta')\n",
    "\n",
    "    Champs = four_point_args(beta0, 0.231+0.1432*1j, (2,3/2,'n') ).channel('s')\n",
    "    \n",
    "    F = ND4(Champs, L = 8, blocklevel = 12) \n",
    "\n",
    "    for x in range(1,5):\n",
    "        s_val = F.values(x/10, channel = True)\n",
    "        t_val = F.values(x/10, channel = False, norm = True)\n",
    "        \n",
    "        print( 'x = %s' % (x/10) )\n",
    "        print( 's-ch: %s' % s_val.real )\n",
    "        print( 't-ch: %s' % t_val.real )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Case of four different fields"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if __name__ == '__main__':\n",
    "    \n",
    "    \"\"\"Functions with four different fields, arbitrary dimensions for diagonal fields.\"\"\"\n",
    "    \n",
    "    c = -0.5432\n",
    "    eps = 0\n",
    "    C = Charge('c', c+eps)\n",
    "    beta0 = C.get('beta')\n",
    "\n",
    "    Champs = four_point_args(beta0, 0.2314, (4,-7/2 ,'n') , 1.265, (2,1/2,'n')).channel('s')\n",
    "    \n",
    "    F = ND4(Champs, L = 12, blocklevel = 20)\n",
    "    \n",
    "    for x in range(1,10):\n",
    "        s_val = F.values(x/10, channel = True)\n",
    "        t_val = F.values(x/10, channel = False, norm = True) # It's possible to turn off t-ch norm.\n",
    "\n",
    "        print( 'x = %s' % (x/10) )\n",
    "        print( 's-ch: %s' % s_val.real )\n",
    "        print( 't-ch: %s' % t_val.real )\n"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}