numerics.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "85c852d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext cython"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94f8aecb",
   "metadata": {},
   "source": [
    "Let $n\\in \\mathbb{N}$, $n\\geq 2$, $p\\in \\mathbb{R}$, $1<p<+\\infty$, $\\mathbb{S}^{n-1}$ be the $n-1$-dimensional unit sphere, and\n",
    "$$\n",
    "  K_{p,n}= |\\mathbb{S}^{n-1}|^{-1} \\int_{\\mathbb{S}^{n-1}}\n",
    "  \\left|e^T s\\right|^{p}\n",
    "  \\,\\mathrm{d}s = \n",
    "  \\frac{1}{\\pi^{1/2}}\\frac{\\Gamma(\\frac{n}{2})\\Gamma(\\frac{p+1}{2})}{\\Gamma(\\frac{n+p}{2})},\n",
    "$$\n",
    "and $w_\\delta$ be a nonnegative radial function supported in $B(0,\\delta)$ normalized in such a way that\n",
    "$$\n",
    "\\int_{B(0,\\delta)} |z|^p w^p_\\delta(z)\\,\\mathrm{d}z = K_{p,n}^{-1}.\n",
    "$$\n",
    "Furthermore, let $\\sigma \\in C^2_c(\\mathbb{R}^n;\\mathbb{R}^n)$. Our goal is to check that\n",
    "$$\n",
    "\\lim_{\\delta\\searrow 0} \\lim_{\\epsilon \\searrow 0}\n",
    "\\int_{B(0,\\delta)\\setminus B(0,\\epsilon)}\n",
    "|\\sigma(x+z)|^{2-p} |\\sigma(x+z)\\cdot z|^{p-2} \n",
    "[\\sigma(x+z)\\cdot z] w^p_\\delta(z)\\,\\mathrm{d}z = \\text{div} \\sigma(x),\n",
    "$$\n",
    "for each $x \\in \\mathbb{R}^n$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b507068",
   "metadata": {},
   "source": [
    "Our main example is going to be $w_\\delta(z) = c_{p,n,\\alpha} |z|^{\\alpha}$, where $c_{p,n,\\alpha}$ is the normalization constant and $-n<p(1+\\alpha)\\leq 0$, i.e., $-1-n/p < \\alpha \\leq -1$. With this in mind we integrate\n",
    "$$\n",
    "\\int_{B(0,\\delta)} |z|^p |z|^{p\\alpha} \\,\\mathrm{d}z\n",
    "= |\\mathbb{S}^{n-1}|\\int_{0}^{\\delta} r^{p(1+\\alpha) + n - 1} \\,\\mathrm{d}r\n",
    "%= |\\mathbb{S}^{n-1}|\\frac{1}{p(1+\\alpha) + n} [r^{p(1+\\alpha) + n}]_{0}^{\\delta} \n",
    "= \\frac{2\\pi^{n/2}}{\\Gamma(\\frac{n}{2})}\\frac{\\delta^{p(1+\\alpha) + n}}{p(1+\\alpha) + n}.\n",
    "$$\n",
    "Thus the normalization constant should be taken as\n",
    "$$\n",
    "c_{p,n,\\alpha} \n",
    "= \\frac{K_{p,n}^{-1}}{\\frac{2\\pi^{n/2}}{\\Gamma(\\frac{n}{2})}\\frac{\\delta^{p(1+\\alpha) + n}}{p(1+\\alpha) + n}}\n",
    "=\\pi^{1/2}\\frac{\\Gamma(\\frac{n+p}{2})}{\\Gamma(\\frac{n}{2})\\Gamma(\\frac{p+1}{2})}\n",
    "\\frac{\\Gamma(\\frac{n}{2})}{2\\pi^{n/2}}\\frac{p(1+\\alpha) + n}{\\delta^{p(1+\\alpha) + n}}\n",
    "=\\frac{\\Gamma(\\frac{n+p}{2})}{2\\pi^{(n-1)/2}\\Gamma(\\frac{p+1}{2})}\n",
    "\\frac{p(1+\\alpha) + n}{\\delta^{p(1+\\alpha) + n}}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d336774",
   "metadata": {},
   "source": [
    "Finally, our main example of $\\sigma$ is going to be a function, whose components are random polynomials of degree $3$. Such functions obviously do not have compact support, but the value of the integral we are trying to compute depends only on the function in the vicinity of $x$, so this is not a problem.  Furthermore, we simply put $x=0$ without any loss of generality."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b768e986",
   "metadata": {},
   "source": [
    "The integrand $$\\sigma_i(x+z) = \\sigma_{0,i} + \\sum_{j} \\sigma_{1,i,j}z_j + \\sum_{j,k} \\sigma_{2,i,j,k}z_jz_k + \\sum_{j,k,\\ell} \\sigma_{3,i,j,k,\\ell}z_jz_kz_\\ell$$ is implemented in Cython to speed up the integration, see below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2cb41fe3",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython -f -c=-O3 -c=-Wno-deprecated-declarations -c=-Wno-unreachable-code-fallthrough\n",
    "#-a\n",
    "\n",
    "from cpython.pycapsule cimport (PyCapsule_New,\n",
    "                                PyCapsule_GetPointer)\n",
    "from cpython.mem cimport PyMem_Malloc,  PyMem_Free\n",
    "from libc.math cimport sin, cos, pow, abs, sqrt, copysign\n",
    "from libc.stdlib cimport malloc, free\n",
    "import scipy\n",
    "from cython cimport boundscheck, wraparound, cdivision, nonecheck\n",
    "\n",
    "@boundscheck(False)\n",
    "@wraparound(False)\n",
    "@nonecheck(False)\n",
    "@cdivision(True)\n",
    "cdef double sph_jac(int n, double* phi_r):\n",
    "    \"\"\"Jacobian of the spherical coordinates\"\"\"\n",
    "    cdef int i\n",
    "    jac = 1.0\n",
    "    for i in range(n-2):\n",
    "        jac *= pow(sin(phi_r[i]),n-2-i)\n",
    "    jac *= pow(phi_r[n-1],n-1)\n",
    "    return jac\n",
    "\n",
    "@boundscheck(False)\n",
    "@wraparound(False)\n",
    "@nonecheck(False)\n",
    "@cdivision(True)\n",
    "cdef void sph2cart(int n, double* phi_r, double* z):\n",
    "    \"\"\"Convert spherical to cartezian coordinates\"\"\"\n",
    "    cdef double sin_prod = 1.0\n",
    "    cdef int i\n",
    "    cdef double r = phi_r[n-1]\n",
    "    for i in range(n-1):\n",
    "        z[i] = r*sin_prod*cos(phi_r[i])\n",
    "        sin_prod *= sin(phi_r[i])\n",
    "    z[n-1] = r*sin_prod\n",
    "\n",
    "@boundscheck(False)\n",
    "@wraparound(False)\n",
    "@nonecheck(False)\n",
    "@cdivision(True)\n",
    "cdef double integrand_main(int n, double* phi_r, void* user_data):\n",
    "    \"\"\"The integrand, written in Cython\"\"\"\n",
    "    # Extract the value of p, power alpha, normalization constant c,\n",
    "    # and polynomial coefficients sig0, sig1, sig2, and sig3\n",
    "    # We also use the end of the data array as temporary buffer for\n",
    "    # computing z and sigma\n",
    "    # Cython uses array access syntax for pointer dereferencing!\n",
    "    cdef double  p    = (<double*>user_data)[0]\n",
    "    cdef double alpha = (<double*>user_data)[1]\n",
    "    cdef double c     = (<double*>user_data)[2]\n",
    "    # coefficients of sigma\n",
    "    cdef double* sig0 = (<double*>user_data)+3\n",
    "    cdef double* sig1 = (<double*>user_data)+3+n\n",
    "    cdef double* sig2 = (<double*>user_data)+3+n+n*n\n",
    "    cdef double* sig3 = (<double*>user_data)+3+n+n*n+n*n*n\n",
    "    # value of sigma\n",
    "    cdef double* sigma = (<double*>user_data)+3+n+n*n+n*n*n+n*n*n*n\n",
    "    cdef double* z     = (<double*>user_data)+3+n+n*n+n*n*n+n*n*n*n+n\n",
    "    #\n",
    "    cdef int i,j,k,l\n",
    "    cdef double r  = phi_r[n-1]\n",
    "    cdef double wd = c*pow(r,p*alpha)\n",
    "    cdef double sigma_abs, sig_dot_z, sig_dot_z_abs, factor, integrand\n",
    "    cdef double* dptr\n",
    "    #\n",
    "    # first, we need to transform from spherical to Cartesian coordinates\n",
    "    sph2cart(n, phi_r, z)\n",
    "    # then we can evaluate the integrand\n",
    "    # we start by calculating sigma\n",
    "    # we start from 3rd order terms to minimize cancellation\n",
    "    dptr = sig3\n",
    "    for i in range(n):\n",
    "        sigma[i] = 0.0\n",
    "        for j in range(n):\n",
    "            for k in range(n):\n",
    "                for l in range(n):\n",
    "                    sigma[i] += dptr[0] * z[j]*z[k]*z[l]\n",
    "                    dptr     += 1\n",
    "    # now 2nd order terms\n",
    "    dptr = sig2\n",
    "    for i in range(n):\n",
    "        for j in range(n):\n",
    "            for k in range(n):\n",
    "                sigma[i] += dptr[0] * z[j]*z[k]\n",
    "                dptr     += 1\n",
    "    # now 1st order terms\n",
    "    dptr = sig1\n",
    "    for i in range(n):\n",
    "        for j in range(n):\n",
    "            sigma[i] += dptr[0] * z[j]\n",
    "            dptr     += 1\n",
    "    # finally constants\n",
    "    for i in range(n):\n",
    "        sigma[i] += sig0[i]\n",
    "    # now we can calculate the length of sigma and dot product with z\n",
    "    sigma_abs = 0.0\n",
    "    sig_dot_z = 0.0\n",
    "    for i in range(n):\n",
    "        sigma_abs += sigma[i]*sigma[i]\n",
    "        sig_dot_z += sigma[i]*z[i]\n",
    "    sig_dot_z_abs = abs(sig_dot_z)\n",
    "    sigma_abs     = sqrt(sigma_abs)\n",
    "    # avoid divisions by zero\n",
    "    if (sigma_abs > 1.0E-12):\n",
    "        # finally, the integrand\n",
    "        factor    = sig_dot_z_abs/sigma_abs\n",
    "        integrand = sigma_abs * copysign(pow(factor,p-1),sig_dot_z) * wd\n",
    "        return integrand * sph_jac(n, phi_r)\n",
    "    else:\n",
    "        return 0.0\n",
    "\n",
    "#\n",
    "# Pack numpy arrays containing the parameters of the integrand\n",
    "#\n",
    "cdef object pack_user_data(int n, double p, double alpha, double c, double[:] sig0, double[:,:] sig1, double[:,:,:] sig2, double[:,:,:,:] sig3):\n",
    "    \"\"\"Wrap data in a PyCapsule for transport.\"\"\"\n",
    "    # Allocate memory\n",
    "    cdef double* data_ptr = <double*> PyMem_Malloc(sizeof(double)*(3+3*n+n*n+n*n*n+n*n*n*n))\n",
    "    data_ptr[0] = p\n",
    "    data_ptr[1] = alpha\n",
    "    data_ptr[2] = c\n",
    "    cdef int i,j,k,l\n",
    "    for i in range(n):\n",
    "        data_ptr[3+i] = sig0[i]\n",
    "    for i in range(n):\n",
    "        for j in range (n):\n",
    "            data_ptr[3+n+i*n+j] = sig1[i,j]\n",
    "    for i in range(n):\n",
    "        for j in range (n):\n",
    "            for k in range(n):\n",
    "                data_ptr[3+n+n*n+i*n*n+j*n+k] = sig2[i,j,k]\n",
    "    for i in range(n):\n",
    "        for j in range (n):\n",
    "            for k in range(n):\n",
    "                for l in range(n):\n",
    "                    data_ptr[3+n+n*n+n*n*n+i*n*n*n+j*n*n+k*n+l] = sig3[i,j,k,l]\n",
    "    # the last 2n doubles is a buffer\n",
    "    return PyCapsule_New(<void*>data_ptr, NULL, free_user_data)\n",
    "\n",
    "cdef void free_user_data(capsule):\n",
    "    \"\"\"Free the memory our value is using up.\"\"\"\n",
    "    PyMem_Free(PyCapsule_GetPointer(capsule, NULL))\n",
    "\n",
    "    \n",
    "def get_low_level_callable_main(int n, double p, double alpha, double c, double[:] sig0, double[:,:] sig1, double[:,:,:] sig2, double[:,:,:,:] sig3):\n",
    "    # scipy.LowLevelCallable expects the function signature to\n",
    "    # appear as the \"name\" of the capsule\n",
    "    func_capsule = PyCapsule_New(<void*>integrand_main,\n",
    "                                 \"double (int, double *, void *)\",\n",
    "                                 NULL)\n",
    "    data_capsule = pack_user_data(n, p, alpha, c, sig0, sig1, sig2, sig3)\n",
    "    \n",
    "    return scipy.LowLevelCallable(func_capsule, data_capsule)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "64040b67",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import integrate\n",
    "from scipy import special"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d41dafc4",
   "metadata": {},
   "source": [
    "This is just a wrap-around the lower level Cython code, which computes the normalization constant before calling the numerical quadrature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b6cebb2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_integral(n,p,alpha,delta,sig0,sig1,sig2,sig3):\n",
    "    ranges = [[0,np.pi]]*(n-2)\n",
    "    ranges.append([0,2*np.pi])\n",
    "    ranges.append([0,delta])\n",
    "    # compute the area of the sphere first\n",
    "    Sn  = 2*pow(np.pi,n/2) / special.gamma(n/2)\n",
    "    # compute Kpn\n",
    "    Kpn = pow(np.pi,-0.5) * special.gamma(n/2) \\\n",
    "          * special.gamma((p+1)/2) / special.gamma((n+p)/2)\n",
    "    # compute normalization constant\n",
    "    pw  = p*(1+alpha)+n\n",
    "    c   = 1/(Kpn * Sn) * pw/pow(delta,pw)\n",
    "    # the main integral\n",
    "    integrand_main = get_low_level_callable_main(n,p,alpha,c,sig0,sig1,sig2,sig3)\n",
    "    I, abserr      = integrate.nquad(integrand_main, ranges)\n",
    "    return I"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6783147a",
   "metadata": {},
   "source": [
    "Finally, we generate four random arrays of coefficients and compare the \"analytical\" value of the divergence with the numerical approximation of the integral for small $\\delta$.  We expect the two to be close to each other."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7800204f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# filter warnings from scipy.integrate.nquad\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "216da86e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n = 3\n",
      "p = 1.500000\n",
      "alpha = -2.666667\n",
      "delta        = 0.000001\n",
      "div sigma(0) = 0.578685\n",
      "integral     = 0.579680\n",
      "error        = 0.000995\n",
      "CPU times: user 6min 8s, sys: 5.67 s, total: 6min 14s\n",
      "Wall time: 6min 15s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "rng = np.random.default_rng()\n",
    "\n",
    "n = 3\n",
    "p = 1.5\n",
    "alpha = (-1-(n-0.5)/p)\n",
    "\n",
    "print(\"n = %d\" % n)\n",
    "assert(n>=2)\n",
    "print(\"p = %f\" % p)\n",
    "assert(p>1)\n",
    "print(\"alpha = %f\" % alpha)\n",
    "assert(p*(1+alpha) <= 0)\n",
    "assert(p*(1+alpha) > -n)\n",
    "\n",
    "# generate random polynomial coefficients\n",
    "# to get sigma(0)=0 put sig0 to zeros\n",
    "#sig0 = rng.normal(size=[n])\n",
    "#sig1 = rng.normal(size=[n,n])\n",
    "#sig2 = rng.normal(size=[n,n,n])\n",
    "#sig3 = rng.normal(size=[n,n,n,n])\n",
    "sig0 = rng.uniform(low=-1.0,high=1.0,size=[n])\n",
    "sig1 = rng.uniform(low=-1.0,high=1.0,size=[n,n])\n",
    "sig2 = rng.uniform(low=-1.0,high=1.0,size=[n,n,n])\n",
    "sig3 = rng.uniform(low=-1.0,high=1.0,size=[n,n,n,n])\n",
    "\n",
    "delta = 1.0E-06\n",
    "print(\"delta        = %f\" % delta)\n",
    "print(\"div sigma(0) = %f\" % np.trace(sig1))\n",
    "I     = compute_integral(n,p,alpha,delta,sig0,sig1,sig2,sig3)\n",
    "err   = np.abs(I - np.trace(sig1))\n",
    "print(\"integral     = %f\" % I)\n",
    "print(\"error        = %f\" % err)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b75b3d24",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}