prec_quadrature.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**prec_quadrature**: Python implementation of the MATLAB quadrature-finding algorithm\n",
    "contained in this directory, so implemented to use high-precision floating point\n",
    "representations\n",
    "\n",
    "The core idea is that overall convergence of the MATLAB algorithm might be limited\n",
    "by basic roundoff error: the problem $\\sum_k w(k) M_j(x_k) = F_j$ has a near-singular Jacobian\n",
    "in the vicinity of known-good solutions, so even small errors in the calculation of\n",
    "the RHS can lead to a stagnated, unconverged iteration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mpmath # Arbitrary-precision floating point\n",
    "from mpmath import mpf # Object for an arbitrary-precision float\n",
    "from mpmath import mp  # Control structure for changing precision\n",
    "import matplotlib # Visualization\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Disable any multithreading inside numpy/MKL\n",
    "import os\n",
    "os.environ['MKL_NUM_THREADS'] = '1'\n",
    "\n",
    "\n",
    "import numpy as np # Arrays for low (machine)-precision floating point values\n",
    "import scipy.spatial\n",
    "import quadfuncs # Functions to generate quadrature moments\n",
    "import calcquad # Functions to evaluate moments and generate the quadrature rule\n",
    "% matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## One-dimensional and Cartesian Product quadrature\n",
    "\n",
    "The core algorithm of the quadrature-finding algorithm is to find a small set of\n",
    "node locations ($\\vec{x}_i$) and weights ($w_i$) where \n",
    "$\\sum_i w_i M_j(\\vec{x}_i) = \\int_\\Omega M_j(\\vec{x}) dA$ for the moment functions $M_j(\\vec{x})$.\n",
    "Computing this requires a sufficiently-accurate right-hand side for each element,\n",
    "and in turn this is most easily given from an impractically large quadrature rule,\n",
    "generated by a Cartesian product."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Moments and orthogonalization\n",
    "\n",
    "Spherical harmonics can be expressed in terms of Cartesian coordinates on the sphere:\n",
    "the zeroth-order harmonic is the constant 1, the first-order harmonics are $(x,y,z)$,\n",
    "the second-order harmonics are $(x^2,xy,y^2,xz,yz)$, and so on. These form the natural \n",
    "basis for a quadrature rule on the sphere, in the same way that polynomials form the \n",
    "natural basis for quadrature on the real line.\n",
    "\n",
    "This has two complications, however:\n",
    "\n",
    "1. There's a difference between a quadrature rule where \"spherical harmonics up to order N\n",
    "are integrated exactly\" and one where \"the error is proportional to $\\Delta x^N$.''\n",
    "2. The spherical harmonics are orthogonal over the entire sphere, but they are most emphatically\n",
    "not orthogonal on any particular element.\n",
    "\n",
    "These complications are related.  Consider an element centered at one of the poles.  There,\n",
    "local variation is described extremely well -- to the proper error estimate -- using only\n",
    "spherical harmonics involving powers $x$ and $y$, so a quadrature rule that integrates only\n",
    "those terms exactly has the right properties.  However, such a quadrature rule will not exactly\n",
    "integrate $z$.\n",
    "\n",
    "In this element, $z = \\sqrt{1-x^2-y^2} \\approx 1-(x^2+y^2)/2$, and it is nearly collinear with\n",
    "the lower-degree spherical harmonics of $x$ and $y$ only.  This causes a serious numerical issue\n",
    "as elements become smaller, and we must resort to orthogonalization and high-precision arithmetic.\n",
    "\n",
    "We accomplish this using a semi-analytical framework.  The ancillary module `quadfuncs` uses SymPy to compute a set of approximation functions $\\phi_n = x^i y^j z - P_n(x,y)$ for a polynomial $P$, such that $\\phi_n$ is orthogonal (over $|x|, |y| \\leq \\epsilon$) to the full set of moment functions that do not contain a 'z' term _and_ $\\phi_k$ for $k < n$.  This process requires the use of high precision for intermediate results, but the result is a set of functions that can calculate $\\phi$ using ordinary doubles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set parameters \n",
    "NTri = 6\n",
    "NSph = 6\n",
    "r0 = 1e-3\n",
    "\n",
    "print('Computing moment functions for order %d x %d' % (NTri,NSph))\n",
    "(funcs_sph, z_poly, z_approx) = quadfuncs.get_sph_funcs(NSph,NTri,r0,debug_funcs=True)\n",
    "funcs_tri = quadfuncs.get_tri_funcs(NTri,r0)\n",
    "print('... complete')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "funcs = funcs_tri + funcs_sph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1d = r0*np.linspace(-1,1,64)\n",
    "y1d = r0*np.linspace(-1,1,64)\n",
    "(x2d,y2d) = np.meshgrid(x1d,y1d,indexing='xy')\n",
    "(moms) = calcquad.eval_moments(x2d,y2d,funcs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NBoxX = np.int(np.ceil(np.sqrt(moms.shape[0])))\n",
    "NBoxY = np.int(np.ceil(moms.shape[0]*1.0/NBoxX))\n",
    "for idx in range(moms.shape[0]):\n",
    "    plt.subplot(NBoxX,NBoxY,idx+1)\n",
    "    plt.pcolormesh(x1d/r0,y1d/r0,moms[idx,:,:].T,\n",
    "                  vmin=-2,vmax=2,\n",
    "                   cmap='RdYlBu_r')\n",
    "    plt.gca().set_xticks([])\n",
    "    plt.gca().set_yticks([])\n",
    "plt.suptitle('Moment funtions');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hexagonal quadrature\n",
    "\n",
    "The above-defined moment functions are now well-conditioned on grids with a length scale of $r_0$.  They're not exactly orthogonal, but the Netwon's Method procedure is capable of dealing with mild amounts of linear dependence.\n",
    "\n",
    "The quadrature-generation algorithm involves taking a candidate set of points and weights and changing it (iteratively) to match the integrals of the moment functions, which are in turn generated by a high-order, tensor-product quadrature rule involving many points.\n",
    "\n",
    "Functions to accomplish this are included in the creatively-named `calcquad` module."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exact quadrature\n",
    "\n",
    "We begin by exploring just what kind of tensor-product quadrature is necessary for sufficient (near numerical-precision) accuray."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mp.dps = 50\n",
    "r0 = mpf('.001')\n",
    "# Create an irregular hexagon, using a LCG to generate \"random\" but deterministic variation\n",
    "xhex = np.array([r0*(mp.cos((2*mp.pi*i)/6) + 0.1*(2*((mp.sqrt(2)*i) % 1) - 1)) \\\n",
    "                     for i in range(6)])\n",
    "yhex = np.array([r0*(mp.sin((2*mp.pi*i)/6) + 0.1*(2*((mp.sqrt(3)*i) % 1) - 1)) \\\n",
    "                     for i in range(6)])\n",
    "\n",
    "# Very high-order quadrature rule\n",
    "(q1x, q1w) = calcquad.gauss(N=50,dps=50)\n",
    "(q2x, q2y, q2w) = calcquad.hex_quad(xhex,yhex,q1x,q1w)\n",
    "q2w = q2w*calcquad.gnom_weight(q2x,q2y)\n",
    "\n",
    "# Evaluate moment functions at these points to get ''exact'' answer\n",
    "moms = calcquad.eval_moments(q2x,q2y,funcs)\n",
    "exact_quad = np.dot(moms,q2w)\n",
    "\n",
    "# Perform a binary search to find the smallest tensor-product quadrature necessary for the\n",
    "# prescribed tolerance\n",
    "low_N = 1\n",
    "high_N = 50\n",
    "try_N = 10\n",
    "current_err = np.inf\n",
    "\n",
    "while (high_N > low_N):\n",
    "    (q1x, q1w) = calcquad.gauss(N=try_N,dps=50)\n",
    "    (q2x, q2y, q2w) = calcquad.hex_quad(xhex,yhex,q1x,q1w)\n",
    "    q2w = q2w*calcquad.gnom_weight(q2x,q2y)\n",
    "    moms = calcquad.eval_moments(q2x,q2y,funcs)\n",
    "    \n",
    "    trial_quad = np.dot(moms,q2w)\n",
    "    current_err = np.max(np.abs(exact_quad-trial_quad)/exact_quad[0])\n",
    "    print('%2d: %e' % (try_N,current_err))\n",
    "    if (current_err > 1e-16): # Guess is too low\n",
    "        low_N = try_N + 1 # Lowest possible value is more than the current test\n",
    "        try_N = (low_N + high_N) // 2\n",
    "    else:\n",
    "        high_N = try_N\n",
    "        try_N = (low_N + high_N) // 2\n",
    "print('Found acceptable quadrature at %d points' % high_N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a 1D quadrature rule of this size\n",
    "(q1x,q1w) = calcquad.gauss(N=10,dps=20)\n",
    "q1x = q1x.astype('double')\n",
    "q1w = q1w.astype('double')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quadrature point rules\n",
    "\n",
    "Next, we want guidelines on how many points are necessary for Nth order quadrature (both with and without exact spherical harmonics) on irregular hexagons.  The absolute lower bound here would be (# functions)/3, since each quadrature point has three degrees of freedom (x,y,weight), but this lower bound is not tight.\n",
    "\n",
    "This test proceeds in two phases.  First, we construct a quadrature rule for the regular spherical hexagon by iterative reduction, beginning from a larger number of candidate quadrature points.  If that process reaches the target (and passes a santity check where quadrature points are not too far outside the hexagon), then we proceed to calculate quadrature rules for 100 randomly-perturbed hexagons.  In some cases, a rule that works for the regular hexagon breaks for perturbed hexagons, perhaps because of a lack of symmetry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r0 = 1e-3\n",
    "quadrules = {} # Dictionary of parameters for acceptable quadrature rules\n",
    "for order in range(6,0,-1):\n",
    "    NTri = order\n",
    "    for NSph in [NTri,0]:\n",
    "        print('Testing quadrature order %d x %d' % (NTri,NSph))\n",
    "        funcs_sph = quadfuncs.get_sph_funcs(NSph,NTri,r0,prec=75)\n",
    "        funcs_tri = quadfuncs.get_tri_funcs(NTri,r0)\n",
    "        funcs = funcs_tri + funcs_sph\n",
    "        \n",
    "        NMin = int(np.ceil(len(funcs)/3.0))\n",
    "        good_rule = False\n",
    "        \n",
    "        for N_targ in [NMin,NMin+1,NMin+2,NMin+3]:\n",
    "            for N_init in [#N_targ,\n",
    "                           N_targ+5,\n",
    "                           len(funcs_sph)+len(funcs_tri)]:\n",
    "                print('   Start (%d), target (%d)' % (N_init, N_targ))\n",
    "                (xreg,yreg,wreg,loop_err,itercount,quad_reg) = calcquad.calc_quad_init(N_init,N_targ,funcs,\n",
    "                                                                                       q1x,q1w,r0)\n",
    "                if (loop_err > 1e-19):\n",
    "                    print('   ... failed (error)')\n",
    "                    continue\n",
    "                if (any(xreg**2+yreg**2 > (1.2*r0)**2)):\n",
    "                    print('   ... failed (sanity)')\n",
    "                    continue\n",
    "                else:\n",
    "                    print('   ... success (%d iter)' % itercount)\n",
    "                # Test adjustment\n",
    "                failed = 0\n",
    "                adjust = 0\n",
    "                initial = 0\n",
    "\n",
    "                iters = 0\n",
    "                for it in range(100):\n",
    "                    xhex = r0*(np.cos(2*np.pi*np.arange(6)/6) + 1e-1*(2*np.random.rand(6) - 1))\n",
    "                    yhex = r0*(np.sin(2*np.pi*np.arange(6)/6) + 1e-1*(2*np.random.rand(6) - 1))\n",
    "\n",
    "                    (xlow,ylow,wlow,loop_err,itercount) = calcquad.calc_hex_quad(xhex,yhex,xreg,yreg,wreg,\n",
    "                                                                                 quad_reg,funcs,q1x,q1w,r0)\n",
    "                    if (loop_err < 1e-14*quad_reg[0]*N_targ**0.5):\n",
    "                        adjust += 1\n",
    "                        iters += itercount\n",
    "                        continue\n",
    "                    loop_err_adj = loop_err/quad_reg[0]\n",
    "                    # In testing, re-initializing the rule for an element seems to never work when ordinary\n",
    "                    # adjustment fails.\n",
    "#                     (xlow,ylow,wlow,loop_err,itercount,_) = calcquad.calc_quad_init(N_init,N_targ,\n",
    "#                                                                                     funcs,q1x,q1w,r0,\n",
    "#                                                                                     xhex,yhex)\n",
    "#                     if (loop_err < 1e-14*quad_reg[0]):\n",
    "#                         initial += 1\n",
    "#                         iters += itercount\n",
    "#                         continue\n",
    "                    loop_err_init = loop_err/quad_reg[0]\n",
    "                    failed += 1\n",
    "                    break\n",
    "                if (failed > 0):\n",
    "                    print('   ... adjustment failed (%d, %e)' % (it, loop_err_adj))\n",
    "                    continue\n",
    "                else:\n",
    "                    print('   ... adjustment success (%d, %d), %d iter' % (adjust, initial,iters))\n",
    "                    good_rule = True\n",
    "                    quadrules[(NTri,NSph)] = (N_init,N_targ)\n",
    "                    \n",
    "#                     plt.figure()\n",
    "#                     plt.plot(np.cos(2*np.pi*np.arange(7)/6),np.sin(2*np.pi*np.arange(7)/6),'k-')\n",
    "#                     plt.plot(np.cos(2*np.pi*np.arange(65)/64),np.sin(2*np.pi*np.arange(65)/64),'b-')\n",
    "#                     plt.plot(xreg/r0,yreg/r0,'r.')\n",
    "#                     plt.show()\n",
    "                                 \n",
    "                    break\n",
    "            if (good_rule): \n",
    "                # Plot moment functions\n",
    "                plt.figure(figsize=(3.5,2.75))\n",
    "                x1d = r0*np.linspace(-1,1,64)\n",
    "                y1d = r0*np.linspace(-1,1,64)\n",
    "                (x2d,y2d) = np.meshgrid(x1d,y1d,indexing='xy')\n",
    "                (moms) = calcquad.eval_moments(x2d,y2d,funcs)\n",
    "                NBoxY = np.int(np.ceil(np.sqrt(moms.shape[0])))\n",
    "                NBoxX = np.int(np.ceil(moms.shape[0]*1.0/NBoxY))\n",
    "                for idx in range(moms.shape[0]):\n",
    "                    ax = plt.subplot(NBoxX,NBoxY,idx+1)\n",
    "                    plt.contourf(x1d/r0,y1d/r0,np.maximum(-2,np.minimum(2,moms[idx,:,:].T)),\n",
    "                                 np.linspace(-2,2,32),\n",
    "                                   cmap='RdYlBu_r')\n",
    "                    plt.gca().set_xticks([])\n",
    "                    plt.gca().set_yticks([])\n",
    "                plt.suptitle('Moment funtions to order %d' % (NTri),fontsize=10);\n",
    "#                 plt.tight_layout()\n",
    "                plt.savefig('figs/moments_%d_%d.pdf' % (NTri,NSph),dpi=150)\n",
    "                plt.show()\n",
    "                break\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "quadrules"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import and process a grid file\n",
    "\n",
    "This section of the notebook imports and processes a text-based grid file.  This grid file opens with a single line containing an integer, then &lt;integer&gt; lines each containing three floating-point values that define a point on the unit sphere.\n",
    "\n",
    "These points are generator locations for the primal grid, which consists of spherical polygons (hexagons and pentagons, typically).  The convex hull of these points defines a dual grid of spherical triangles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def build_dual_centers(chull,primal_generators):\n",
    "    '''build_dual_centers: calculate the centers of a dual grid (triangles) given primal generator \n",
    "    locations and a convex hull'''\n",
    "\n",
    "    # Build barycenters and voronoi centers of the dual triangles\n",
    "\n",
    "    dual_bcent = np.zeros((len(chull.simplices),3),dtype=np.double)\n",
    "    dual_vcent = np.zeros((len(chull.simplices),3),dtype=np.double)\n",
    "\n",
    "    for (idx, tri) in enumerate(dual_faces):\n",
    "        triverts = primal_generators[tri,:]\n",
    "\n",
    "        # The barycenter is the average of the triangle vertices, projected to the sphere\n",
    "        bcent = np.mean(triverts,axis=0)\n",
    "        bcent = bcent/np.sum(bcent**2)**0.5\n",
    "        dual_bcent[idx,:] = bcent\n",
    "\n",
    "        # The voronoi center is the point orthogonal to the triangle edges; see Muira [2005], equation 2\n",
    "        vcent = np.cross(triverts[1,:]-triverts[0,:],triverts[2,:]-triverts[0,:])\n",
    "        vcent = vcent/np.sum(vcent**2)**0.5 * np.sign(np.dot(vcent,triverts[0,:]))\n",
    "        dual_vcent[idx,:] = vcent\n",
    "        \n",
    "    return(dual_bcent,dual_vcent)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Build list of vertices (by face) for the primal grid\n",
    "def build_verts(primal_generators,dual_faces,dual_centers):\n",
    "    '''build_verts: Build a list of vertices (triangle centers) for the primal grid (hexagon/pentagon)'''\n",
    "    primal_vertices = [ [] for x in range(len(primal_generators))]\n",
    "\n",
    "    # Iterate through dual grid vertices, associating each face with its vertex (primal face)\n",
    "    for (idx, tri) in enumerate(dual_faces):\n",
    "        for t in tri:\n",
    "            primal_vertices[t].append(idx)\n",
    "\n",
    "    # Now, go through this list and order the vertices in CCW direction\n",
    "    for (idx, verts) in enumerate(primal_vertices):\n",
    "        vert_coords = dual_centers[verts,:]\n",
    "\n",
    "        # Compute an approximate center of the polygon\n",
    "        poly_center = np.mean(vert_coords,axis=0)\n",
    "\n",
    "        # Generate false east/north coordinates, to determine \n",
    "        false_east = vert_coords[0,:] - poly_center\n",
    "        false_north = np.cross(poly_center,false_east)\n",
    "\n",
    "        false_x = np.dot(false_east,vert_coords.T)\n",
    "        false_y = np.dot(false_north,vert_coords.T)\n",
    "\n",
    "        angle = np.mod(np.arctan2(false_y,false_x)+2*np.pi,2*np.pi)\n",
    "        angle[0] = 0.0\n",
    "\n",
    "        sorted_idx = np.argsort(angle)\n",
    "        primal_vertices[idx] = [verts[i] for i in sorted_idx]\n",
    "    return primal_vertices\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now enumerate the edges of this grid.  Each primal edge (edge of a spherical polygon) has a dual edge (edge of a spherical triangle), so it is defined uniquely by the 4-tuple (primal face 1, primal face 2, dual face 1, dual face 2).  From the data structure given by the ConvexHull, we can build this tuple in two passes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def build_edges(dual_faces):\n",
    "    '''build_edges: Build a list of edges (face/vertex connections)'''\n",
    "    edgedict = {} # Dictionary indexed by (primal face 1, primal face 2)\n",
    "    def addvert(edge, pf1, pf2, df1):\n",
    "        # For the given dual edge connecting (pf1) and (pf2), add \n",
    "        # the connection (df1) to the edge dictionary\n",
    "\n",
    "        val = edge.get((pf1,pf2),[])\n",
    "        val.append(df1)\n",
    "        edge[(pf1,pf2)]=sorted(val)\n",
    "\n",
    "    # Pass 1 -- build edges by accumulating primal face / dual face associations\n",
    "    for (idx,tri) in enumerate(dual_faces):\n",
    "        stri = sorted(tri) # Sort vertices of the triangle to give a canonical order\n",
    "        addvert(edgedict,stri[0],stri[1],idx)\n",
    "        addvert(edgedict,stri[1],stri[2],idx)\n",
    "        addvert(edgedict,stri[0],stri[2],idx)\n",
    "\n",
    "    # Pass 2 -- build an edge array, where each entry contains:\n",
    "    # (pf1, pf2, pv1, pv2)\n",
    "    # to represent the edge between faces pf1/pf2 that consists of vertices pv1/pv2\n",
    "    edges = np.zeros((len(edgedict),4),dtype=np.int)\n",
    "    for (idx,(k, v)) in enumerate(edgedict.iteritems()):\n",
    "        if (v[0] < v[1]):\n",
    "            edges[idx,:] = [k[0], k[1], v[0], v[1]]\n",
    "        else:\n",
    "            edges[idx,:] = [k[0], k[1], v[1], v[0]]\n",
    "\n",
    "    # Sort the edges array so that faces are in lexical order\n",
    "    edges = edges[edges[:,1].argsort()]\n",
    "    edges = edges[edges[:,0].argsort()]\n",
    "    return edges"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define test functions to evaluate the accuracy of the resulting quadrature rules.  `test_sph0` through `test_sph3` are polynomials in Cartesian coordinates, which are also linear combinations of spherical harmonics.  Quadrature rules which include the \"spherical\" moments (of sufficiently high order) should integrate these functions exactly, up to the prescribed tolerance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_tanh(xyz):\n",
    "    '''test_tanh: tanh test function (Reeger and Fornberg 2015, eqn f_1)'''\n",
    "    import numpy as np\n",
    "    return (1+np.tanh(9*(xyz[:,2]-xyz[:,0]-xyz[:,1])))\n",
    "\n",
    "def test_sph0(xyz):\n",
    "    '''test_sph0 -- cell area'''\n",
    "    return(xyz[:,0]**0)\n",
    "\n",
    "def test_sph1(xyz):\n",
    "    '''test_sph1: first-order spherical harmonics'''\n",
    "    return (xyz[:,0]+xyz[:,2])\n",
    "\n",
    "def test_sph2(xyz):\n",
    "    '''test_sph2: second-order spherical harmonics'''\n",
    "    return (xyz[:,0]**2 + xyz[:,0]*xyz[:,2])\n",
    "\n",
    "def test_sph3(xyz):\n",
    "    '''test_sph3: third-order spherical harmonics'''\n",
    "    return xyz[:,0]**3 + xyz[:,1]**3 + xyz[:,0]*xyz[:,1]*xyz[:,2]\n",
    "\n",
    "test_funcs = ['1','$x+z$','$x^2 + xz$','$x^3 + y^3 + xyz$',\n",
    "    '$1+\\\\tanh(9(z-x-y))$']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need a sufficiently accurate quadrature rule to generate reference values for the test integrals, so copy the tensor-product generation code from `calcquad` and apply the rule directly to provided functions.  There is no need to save the resulting quadrature rule, and in fact doing so would chew through memory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def control_quad(verts,q1x,q1w,funcs):\n",
    "    '''control_quad: calculate tensor-product quadrature of provided functions over a single cell'''\n",
    "    \n",
    "    import numpy as np\n",
    "    import calcquad\n",
    "    \n",
    "    # Calculate an approximate center of the cell, by averaging the vertices\n",
    "    center = np.mean(verts,axis=0)\n",
    "    # Project back to the unit sphere\n",
    "    center /= np.linalg.norm(center)\n",
    "\n",
    "    if (verts.shape[0] == 5):\n",
    "        # We have a hexagon, not a pentagon, so duplicate a vertex\n",
    "        verts = verts[[0]+range(5),:]\n",
    "\n",
    "    t_vec = center # t is the outward coordinate, which takes the place of 'z'\n",
    "    # Pick one of the vertex directions as 'r' (x-substitute)\n",
    "    r_vec = verts[0,:] - center\n",
    "    r_vec -= t_vec*np.dot(r_vec,t_vec)\n",
    "    r_vec /= np.linalg.norm(r_vec)\n",
    "\n",
    "    # Pick the other direction as 's' (y-substitute)\n",
    "    s_vec = np.cross(t_vec,r_vec)\n",
    "\n",
    "    # These vectors give a transform matrix\n",
    "    tmat = np.array([r_vec,s_vec,t_vec]) # Convert x/y/z to r/s/t\n",
    "    tmati = np.linalg.inv(tmat) # Convert r/s/t to x/y/z\n",
    "\n",
    "    # Convert the hexagon vertices to {r,s,t} space\n",
    "    proj_verts = np.dot(tmat,verts.T)\n",
    "\n",
    "    # Divide the coordinates by their third component to effect the gnomonic projection\n",
    "    hex_gr = proj_verts[0,:]/proj_verts[2,:]\n",
    "    hex_gs = proj_verts[1,:]/proj_verts[2,:]\n",
    "\n",
    "    (quad_gr,quad_gs,quad_w) = calcquad.hex_quad(hex_gr,hex_gs,q1x,q1w)\n",
    "    quad_w *= calcquad.gnom_weight(quad_gr,quad_gs)\n",
    "    \n",
    "    # Convert the quadrature points on the gnonomic (r,s) plane back to a rotated sphere:\n",
    "\n",
    "    # First, place the points on the tangent plane\n",
    "    quad_verts = np.dot(tmat.T,np.vstack((quad_gr, quad_gs, quad_gr**0)))\n",
    "\n",
    "    # Then, normalize to the unit sphere\n",
    "    quad_verts = quad_verts / np.sqrt(np.sum(quad_verts**2,axis=0))\n",
    "    quad_verts = quad_verts.T\n",
    "\n",
    "    \n",
    "    out = map(lambda f: np.dot(quad_w,f(quad_verts)), funcs)\n",
    "    return out\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_ortho(verts,points,center,east):\n",
    "    '''Helper function to plot a polygon and points in an orthographic projection'''\n",
    "    \n",
    "    # Vectors tangent to the sphere\n",
    "    t_vec = center\n",
    "    t_vec /= np.linalg.norm(t_vec)\n",
    "    r_vec = east - np.dot(east,t_vec)*t_vec\n",
    "    r_vec /= np.linalg.norm(r_vec)\n",
    "    s_vec = np.cross(t_vec,r_vec)\n",
    "    tmat = np.array([r_vec,s_vec,t_vec])\n",
    "    \n",
    "    for idx in range(-1,verts.shape[0]-1):\n",
    "        v1 = verts[idx,:]\n",
    "        v2 = verts[idx+1,:]\n",
    "        \n",
    "        gnom_1 = np.dot(tmat,v1)\n",
    "        gnom_1 /= gnom_1[2]\n",
    "        gnom_2 = np.dot(tmat,v2)\n",
    "        gnom_2 /= gnom_2[2]\n",
    "        \n",
    "        ignom = gnom_1[:,np.newaxis] + (gnom_2-gnom_1)[:,np.newaxis]*np.linspace(0,1,32)\n",
    "        iortho = ignom / (1+ignom[0,:]**2 + ignom[1,:]**2)**0.5\n",
    "        \n",
    "        plt.plot(iortho[0,:],iortho[1,:],'k-')\n",
    "        \n",
    "    proj_pts = np.dot(tmat,points.T)\n",
    "    plt.plot(proj_pts[0,:],proj_pts[1,:],'r.')\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Scaffolding for parallel generation of reference results\n",
    "import ipyparallel\n",
    "try:\n",
    "    rc.queue_status()\n",
    "except (NameError, IOError):\n",
    "    rc = ipyparallel.Client()\n",
    "    view = rc[:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Loop over all the admissible quadrature rules and over a selection of grid refinement \n",
    "# levels to produce quadrature rules and error estimates\n",
    "qtensor_dict = {}\n",
    "qcalc_dict = {}\n",
    "verts_dict = {}\n",
    "for (NTri,NSph) in sorted(quadrules.keys()):\n",
    "    quad_tensor = []\n",
    "    quad_calc = []\n",
    "    qp_arr = []\n",
    "    qw_arr = []\n",
    "    (N_start, N_targ) = quadrules[(NTri,NSph)]\n",
    "    for LEVEL in range(1,7):\n",
    "        gridfile = '../iModel/grid/icos_pol_scvt_h1_%d.xyz' % LEVEL\n",
    "\n",
    "        with open(gridfile,'r') as genfile:\n",
    "            npoints = int(genfile.readline())\n",
    "            primal_generators = np.array([[float(f) for f in l.split()] for l in genfile])\n",
    "\n",
    "        chull = scipy.spatial.ConvexHull(primal_generators)\n",
    "        dual_faces = chull.simplices # Triangles, indexed by generator\n",
    "        (dual_bcent,dual_vcent) = build_dual_centers(chull,primal_generators)\n",
    "        primal_vertices = build_verts(primal_generators,dual_faces,dual_vcent)\n",
    "        edges = build_edges(dual_faces)\n",
    "\n",
    "        (quad_points, quad_weights, rel_errs, total_iter,\n",
    "         xreg, yreg, wreg) = calcquad.quadgrid(primal_vertices,dual_vcent,q1x,q1w,NTri,NSph,N_start,N_targ)\n",
    "\n",
    "        qp_arr.append(np.array(quad_points))\n",
    "        qw_arr.append(np.array(quad_weights))\n",
    "\n",
    "        # Get control quadratures\n",
    "        verts = [dual_vcent[primal_vertices[i],:] for i in range(len(primal_vertices))]\n",
    "        if LEVEL not in verts_dict:\n",
    "            verts_dict[LEVEL] = verts\n",
    "        view.push(dict(control_quad = control_quad,\n",
    "                       test_tanh=test_tanh,\n",
    "                       test_sph0 = test_sph0,\n",
    "                       test_sph1 = test_sph1,\n",
    "                       test_sph2 = test_sph2,\n",
    "                       test_sph3 = test_sph3,\n",
    "                       q1x=q1x,q1w=q1w))\n",
    "        quad_tensor.append(np.array(view.map_sync(lambda v: \n",
    "                                                  control_quad(v,q1x,q1w,\n",
    "                                                               [test_sph0,test_sph1,test_sph2,test_sph3,test_tanh]), \n",
    "                                                  verts)))\n",
    "\n",
    "        qp2d = np.array(quad_points).reshape((-1,3))\n",
    "        qw2d = np.array(quad_weights)\n",
    "\n",
    "        quadcalc_arr = np.array(\n",
    "            map(lambda f: np.sum(f(qp2d).reshape((-1,N_targ))*qw2d,axis=1),\n",
    "                [test_sph0,test_sph1,test_sph2,test_sph3,test_tanh])).T\n",
    "        quad_calc.append(quadcalc_arr)\n",
    "\n",
    "        if (LEVEL == 2):\n",
    "            # Plot initial and final quadrature rules at a sample hexagon\n",
    "            r0 = (4.0/len(primal_vertices))**0.5 \n",
    "            xhex_g = r0*np.cos(np.arange(6)*2*np.pi/6) # Regular hexagon\n",
    "            zhex_g = r0*np.sin(np.arange(6)*2*np.pi/6)\n",
    "            yhex = (1+xhex_g**2 + zhex_g**2)**-0.5\n",
    "            vhex = np.array([xhex_g*yhex,yhex,zhex_g*yhex]).T\n",
    "            # Regular quadrature points, projected to sphere\n",
    "            zreg = (1+xreg**2+yreg**2)**-0.5\n",
    "            ax = plt.figure(figsize=(3.5,2.75))\n",
    "            plot_ortho(vhex,np.array((xreg*zreg,zreg,yreg*zreg)).T,np.array([0.0,1,0]),np.array([1.0,0,0]))\n",
    "            plt.title('Regular hexagon (%d, %d)' % (NTri,NSph))\n",
    "            plt.savefig('figs/regquad_%d_%d.pdf' % (NTri,NSph), dpi=150)\n",
    "            plt.show()\n",
    "            IDX = 12\n",
    "            ax = plt.figure(figsize=(3.5,2.75))\n",
    "            cent = np.mean(verts[IDX],axis=0)\n",
    "            plot_ortho(verts[IDX],quad_points[IDX],cent,verts[IDX][0,:] - cent)\n",
    "            plt.title('Sample cell (%d, %d)' % (NTri,NSph))\n",
    "            plt.savefig('figs/sampquad_%d_%d.pdf' % (NTri,NSph), dpi=150)\n",
    "            plt.show()\n",
    "\n",
    "\n",
    "        # Write out 0-indexed grid information\n",
    "\n",
    "        if (False):\n",
    "            # Write vertex data: each entry is the (x,y,z) coordinate of a primal-grid vertex\n",
    "            with open('vertices_%d.xyz' % LEVEL,'w') as vertfile:\n",
    "                # One header line: number of vertices\n",
    "                vertfile.write('%d\\n' % dual_vcent.shape[0])\n",
    "                for idx in range(dual_vcent.shape[0]):\n",
    "                    vertfile.write(' '.join(['%+.17e' % dual_vcent[idx,j] for j in range(3)]) + '\\n')\n",
    "\n",
    "\n",
    "            # Write edge data: each entry is the (v1, v2, f1, f2) tuple (zero-indexed) of\n",
    "            # faces and vertices corresponding to a grid edge\n",
    "            with open('edges_%d.ij' % LEVEL,'w') as edgefile:\n",
    "                edgefile.write('%d\\n' % edges.shape[0])\n",
    "                for idx in range(edges.shape[0]):\n",
    "                    edgefile.write('%d %d %d %d\\n' % tuple(edges[idx,:]))\n",
    "\n",
    "            # Write face data: each entry is the (v0..v[45]) tuple (zero-indexed) of vertices\n",
    "            # corresponding to a grid face, in anticlockwise order\n",
    "            with open('faces_%d.ij' % LEVEL,'w') as facefile:\n",
    "                facefile.write('%d\\n' % len(primal_vertices))\n",
    "                for vlist in primal_vertices:\n",
    "                    facefile.write(' '.join('%d' % v for v in vlist) + '\\n')\n",
    "\n",
    "            with open('quad_%d_%dx%d.xyzw' % (LEVEL,NTri,NSph),'w') as qfile:\n",
    "                qfile.write('%d %d\\n' % (len(quad_weights), len(quad_weights[0])))\n",
    "                for (idx) in range(len(quad_points)):\n",
    "                    for (qp,qw) in zip(quad_points[idx],quad_weights[idx]):\n",
    "                        qfile.write('%+.17e %+.17e %+.17e %+.17e\\n' % (tuple(qp)+(qw,)))\n",
    "                        \n",
    "    qtensor_dict[(NTri,NSph)] = quad_tensor\n",
    "    qcalc_dict[(NTri,NSph)] = quad_calc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot error convergence curves\n",
    "\n",
    "lines = []\n",
    "labels = []\n",
    "plt.figure(figsize=(3.5,2.75))\n",
    "for NTri in range(1,7):\n",
    "    NSph = 0\n",
    "    quad_calc = qcalc_dict[(NTri,NSph)]\n",
    "    quad_tensor = qtensor_dict[(NTri,NSph)]\n",
    "    (N_start,N_targ) = quadrules[(NTri,NSph)]\n",
    "    diffs = []\n",
    "    for idx in range(len(quad_calc)):\n",
    "        quad_diff = (quad_calc[idx] - quad_tensor[idx]).T/quad_calc[idx][:,0]\n",
    "        diffs.append(np.max(np.abs(quad_diff),axis=1))\n",
    "#         print(diffs[-1])\n",
    "    diffs = np.array(diffs)\n",
    "    pts = np.array([g.shape[0] for g in quad_calc],dtype=np.double)*N_targ\n",
    "#     plt.loglog(pts,diffs[:,1],linestyle=':',marker='.',color=plt.cm.tab10(NTri-1))\n",
    "#     plt.loglog(pts,diffs[:,2],linestyle='--',marker='.',color=plt.cm.tab10(NTri-1))\n",
    "#     plt.loglog(pts,diffs[:,3],linestyle='-.',marker='.',color=plt.cm.tab10(NTri-1))\n",
    "    ln, = plt.loglog(pts,diffs[:,4],linestyle='-',marker='.',color=plt.cm.tab10(NTri-1))\n",
    "    lines.append(ln)\n",
    "    labels.append('Order %d' % (NTri))\n",
    "\n",
    "    \n",
    "plt.xlabel('# Points')\n",
    "plt.ylabel('$||\\cdot||_\\infty$ error')\n",
    "plt.legend(lines,labels)\n",
    "plt.title('Quadrature errors (triangular)')\n",
    "plt.axis((42,800000,3e-10,1e0))\n",
    "plt.tight_layout()\n",
    "plt.savefig('figs/tri_err.pdf',dpi=150)\n",
    "plt.show()\n",
    "\n",
    "lines = []\n",
    "labels = []\n",
    "plt.figure(figsize=(3.5,2.75))\n",
    "for NTri in range(1,6):\n",
    "    NSph = NTri\n",
    "    quad_calc = qcalc_dict[(NTri,NSph)]\n",
    "    quad_tensor = qtensor_dict[(NTri,NSph)]\n",
    "    (N_start,N_targ) = quadrules[(NTri,NSph)]\n",
    "    diffs = []\n",
    "    for idx in range(len(quad_calc)):\n",
    "        quad_diff = (quad_calc[idx] - quad_tensor[idx]).T/quad_calc[idx][:,0]\n",
    "        diffs.append(np.max(np.abs(quad_diff),axis=1))\n",
    "#         print(diffs[-1])\n",
    "    diffs = np.array(diffs)\n",
    "    pts = np.array([g.shape[0] for g in quad_calc],dtype=np.double)*N_targ\n",
    "#     plt.loglog(pts,diffs[:,1],linestyle=':',marker='.',color=plt.cm.tab10(NTri-1))\n",
    "#     plt.loglog(pts,diffs[:,2],linestyle='--',marker='.',color=plt.cm.tab10(NTri-1))\n",
    "#     plt.loglog(pts,diffs[:,3],linestyle='-.',marker='.',color=plt.cm.tab10(NTri-1))\n",
    "    ln, = plt.loglog(pts,diffs[:,4],linestyle='-',marker='.',color=plt.cm.tab10(NTri-1))\n",
    "    lines.append(ln)\n",
    "    labels.append('Order %d' % (NTri))\n",
    "\n",
    "    \n",
    "plt.xlabel('# Points')\n",
    "plt.ylabel('$||\\cdot||_\\infty$ error')\n",
    "plt.legend(lines,labels)\n",
    "plt.title('Quadrature errors (complete)')\n",
    "plt.axis((42,800000,3e-10,1e0))\n",
    "plt.tight_layout()\n",
    "plt.savefig('figs/full_err.pdf',dpi=150)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot a sample visualization of quadrature error.  For the pre-selected values and the tanh\n",
    "# test function, the quadrature error should be concentrated in the small transition band\n",
    "\n",
    "polys = []\n",
    "colors = []\n",
    "pts = []\n",
    "NTri = 4; NSph = 4\n",
    "SIZE_IDX = -1\n",
    "quad_calc = qcalc_dict[(NTri,NSph)]\n",
    "quad_tensor = qtensor_dict[(NTri,NSph)]\n",
    "quad_diff = ((quad_calc[SIZE_IDX] - quad_tensor[SIZE_IDX]).T/quad_calc[SIZE_IDX][:,0]).T\n",
    "cmap = plt.get_cmap('RdYlBu')\n",
    "TESTFUNC=4\n",
    "\n",
    "# Plot true function\n",
    "vmax = np.max(np.abs(quad_calc[SIZE_IDX][:,TESTFUNC]/quad_calc[SIZE_IDX][:,0]))\n",
    "from matplotlib.collections import PatchCollection\n",
    "from matplotlib.patches import Polygon\n",
    "for (idx,v) in enumerate(verts):\n",
    "    if (np.any(v[:,1] < 0)): continue\n",
    "    polys.append(Polygon(v[:,[0,2]],closed=True))\n",
    "    colors.append(quad_calc[SIZE_IDX][idx,TESTFUNC]/quad_calc[SIZE_IDX][idx,0])\n",
    "fig, ax = plt.subplots()\n",
    "p = PatchCollection(polys, cmap=plt.get_cmap('RdYlBu_r'), alpha=1.0)\n",
    "p.set_array(np.array(colors))\n",
    "p.set_clim([-vmax,vmax])\n",
    "p.set_edgecolor(None)\n",
    "plt.colorbar(p)\n",
    "ax.add_collection(p)\n",
    "ax.autoscale()\n",
    "ax.set_title(test_funcs[TESTFUNC])\n",
    "\n",
    "# pos_pts = (qp2d[:,1] > 0)\n",
    "# plt.plot(qp2d[pos_pts,0],qp2d[pos_pts,2],'k.',markersize=4)\n",
    "plt.show()\n",
    "\n",
    "# Plot quadrature difference\n",
    "vmax = np.max(np.abs(quad_diff[:,TESTFUNC]))\n",
    "cmap = plt.get_cmap('RdYlBu')\n",
    "from matplotlib.collections import PatchCollection\n",
    "from matplotlib.patches import Polygon\n",
    "for (idx,v) in enumerate(verts):\n",
    "    if (np.any(v[:,1] < 0)): continue\n",
    "    polys.append(Polygon(v[:,[0,2]],closed=True))\n",
    "    colors.append(quad_diff[idx,TESTFUNC])\n",
    "    pts.append(quad_points[idx])\n",
    "fig, ax = plt.subplots()\n",
    "p = PatchCollection(polys, cmap=plt.get_cmap('RdYlBu_r'), alpha=1.0)\n",
    "p.set_array(np.array(colors))\n",
    "p.set_clim([-vmax,vmax])\n",
    "p.set_edgecolor(None)\n",
    "plt.colorbar(p)\n",
    "ax.add_collection(p)\n",
    "ax.autoscale()\n",
    "ax.set_title('%s error' % test_funcs[TESTFUNC])\n",
    "pts = np.array(pts).reshape((-1,3))\n",
    "# plt.plot(pts[:,0],pts[:,2],'k.',markersize=2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}