clean up and notes

src/madness/bspline/bspline.cc

@@ -613,8 +613,8 @@ class BsplineBasis : protected knotsT {
 
     // Test the basic class functionality --- this is not yet a unit test
     static bool test() {
-        size_t nknots = 50;
-        size_t order = 12; 
+        size_t nknots = 15;
+        size_t order = 6; 
         T xlo = 0;//1e-5; //-1; //-constants::pi;
         T xhi = 25;  //2*constants::pi;
         T xtest = 1; //0.5*(xhi+xlo);
@@ -637,7 +637,12 @@ class BsplineBasis : protected knotsT {
         // }
 
         //print(B.tabulate_basis(knots));
-        //print(B.tabulate_basis(xsam));
+        {
+            auto xxsam = oversample_knots(oversample_knots(oversample_knots(oversample_knots(xsam))));
+            auto A = B.tabulate_basis(xxsam);
+            Gnuplot G("set style data lines; set title 'B-spline basis'; set xlabel 'x'; set ylabel 'b[i](x)'; set key off","basis.gnuplot");
+            gnuplot_tensor(G, xxsam, A);
+        }
 
         Tensor<T> f(nsam), df(nsam), d2f(nsam);
         for (size_t i=0; i<nsam; i++) {
@@ -1261,6 +1266,23 @@ class BsplineFunction {
     }
 };
 
+// The spherical component of the GF is 2 (H(r-s) h(r) j(s) + H(s-r) h(s) j(r)) with
+// * H the Heaviside function
+// * h(r) = exp(-mu*r)/r
+// * j(s) = sinh(mu*s)/(mu*s) = (exp(mu*s) - exp(-mu*s))/(2*mu*s)
+// When applying to a function f(r) the first term gives the integral
+// (1/mu r) int_0^r s (exp(-mu*(r-s)) - exp(-mu*(r+s))) f(s) ds
+// and the second term gives the integral
+// (1/mu r) int_r^inf s (exp(-mu*(s-r)) - exp(-mu*(s+r))) f(s) ds
+// In both of the these the arguments to the exponentials are negative so nothing blows up.
+// These are again computed using GL quadrature.
+// int_0^r g(r) = \sum^{whole knot intervals} + \sum^{partial knot interval}
+// whole knot intervals are easy because they reuse the current quadrature points
+// \int_{interval} = sum_{quadrature points} w_i g(x_i) with w_i and x_ from XW as usual
+// For a partial interval, we need additional quadrature points since the range is now
+// [lo,r] or [r,hi].  Note you dont need both since [r,hi] = [lo,hi] - [lo,r], the former being the full interval.
+// 
+
 
 // // Make the matrix that applies the second derivative after projecting into a a two higher order basis.
 // // Gives us 1 order higher accuracy and produces a result in the same order basis as the input.
@@ -1321,7 +1343,7 @@ int main() {
     //BsplineBasis<double,KnotsUniform<double>>::test();
     //BsplineBasis<double,KnotsChebyshev<double>>::test();
     //BsplineBasis<double,KnotsGeometricShifted<double>>::test();
-    //BsplineBasis<double,KnotsRational<double>>::test();
+    BsplineBasis<double,KnotsRational<double>>::test();
     BsplineFunction<double>::test();
     return 0;
 }

src/madness/bspline/notes.tex

@@ -0,0 +1,280 @@
+\documentclass[12pt]{article}
+\usepackage{cite}
+\usepackage{mathtools} %% includes amsmath
+\usepackage{amssymb,amsfonts}
+\usepackage{algorithmic}
+\usepackage{textcomp}
+\usepackage{xcolor}
+\usepackage{physics}
+\usepackage{graphicx,graphics}
+\setlength{\oddsidemargin}{0.25in}
+\setlength{\evensidemargin}{-0.25in}
+\setlength{\topmargin}{-.5in}
+\setlength{\textheight}{9in}
+\setlength{\parskip}{.1in}
+\setlength{\parindent}{2em}
+\setlength{\textwidth}{6.25in}
+
+\newcommand{\N}{\mathcal{N}}
+
+\title{Notes on b-spline implementation}
+\date{\today}
+\author{Robert J. Harrison}
+
+\begin{document}
+
+\maketitle
+
+This document is initially looking at the b-spline only implementation. 
+
+The functional form is
+\begin{eqnarray}
+  f(\vb{r}) & = & \sum_{l m} \N_{lm}(\hat{\vb{r}}) R_{lm}(r) \\
+  R_{lm}(r) & = & \sum_i c_{lmi} b_i(r)
+\end{eqnarray}
+Comments:
+\begin{itemize}
+\item The usual spherical coordinate system with $r = |\vb(r)|$, $\theta$ ($z=r \cos \theta$), $\phi$.
+\item The $\N_{lm}(\hat{\vb{r}})$ where $\hat{r}$ is the unit vector are the real solid spherical harmonics of Yang normalized so that
+\begin{eqnarray}
+  \int_0^{2 \pi} d\phi \int_0^\pi d\theta \sin \theta  \N_{lm} (\hat{r}) \N_{l^\prime m^\prime}(\hat{r})  & = & \delta_{l l^\prime} \delta_{m m^\prime}
+\end{eqnarray}
+Note that Yang does not normalize his harmonics for ease of computing various quantities.  However, we need radial and angular components to be normalized so we can perform numerical thresholding without having to keep track of normalization constants, which for Yang's harmonics (denoted $N_{lm}$)  can be truly huge ($O(10^{2l})$ or worse).  Yang's normalization also results in $N_{10}=z$ but $N_{11}=-x/2$ and $N_{1-1}=-y/2$ and so potentially breaks symmetry.  
+\item $b_i(r)$ is a b-spline basis function defined on knots that include the origin. 
+\item For non-zero angular momentum $l$, we expect the radial functions to be of the form $r^l (a + b r + \cdots)$, thus to exactly represent these functions near the origin the degree of the b-spline should be $p \ge l$.
+\end{itemize}
+
+\section{B-splines}
+
+The basis is specified by 
+\begin{itemize}
+\item $p$ --- the degree of the basis, and
+\item a vector of unique knots.
+\end{itemize}
+We employ the full set of b-splines over a vector of unique knots.  The definition of the full basis requires a vector with the knots padded on either end by the last knot repeated $p$ times.  Note we actually store this vector ($t$ in the code) since it makes algorithms to compute with b-splines much easier. 
+
+Thus, for a basis of degree $p$ (order $p+1$) and $m$ unique knots 
+\begin{itemize}
+\item they form a complete basis over the interval with approximation error $O\left(h^{r+p+1}\right)$ for a function with $r$ continuous derivatives 
+\item the basis has $C_{p-1}$ continuity, meaning its value and all derivatives up to $p-1$ are continuous
+\item the number of padded knots is $m + 2p$
+\item the total number of basis functions is $m + p - 1$ 
+\item the number of basis functions non-zero over a given knot interval is $p + 1$
+\item each spline is non-zero over $p+1$ knot intervals and is controlled by $p+2$ knot values
+\item each spline is positive within its support
+\item the sum of all splines at each knot is unity by construction
+\item the fully interior sides of splines are zero at the end of their support along with their derivatives up to $p-1$
+\item splines touching the edge have lower number of vanishing derivatives, and only the first/last splines have non-zero values at the end points (and these values will be one by construction).
+\end{itemize}
+
+\section{Radial knots}
+
+We are computing on the interval $[0,L]$ which we can obtain by scaling the unit interval $[0,1]$ by $L$.  By defining monotonicaly-increasing maps from $[0,1]$ back onto $[0,1]$ we can nest transformations.    For $n$ knots
+\begin{itemize}
+\item Uniform knots $x(i) = i/(n-1)$
+\item Chebyshev-like knots (that include the end points and concentrate quadratically on both ends) $x(i) = (1-\cos \pi i / (n-1))/2$
+\item Half Chebyshev-like knots (that only concentrate on the left) $x(i) = 1 - \cos \pi i / (2(n-1))$
+\item Rational function similar to half-Cheby $x(i) = i^2 (1+a) / (1 + ai)$
+\item Geometric
+\item Power
+\end{itemize}
+
+The below figures display the basis with 15 knots, order 6, with uniform and rational knots.
+
+\includegraphics[width=.9\linewidth]{basis.pdf}
+
+\includegraphics[width=.9\linewidth]{basis-non-uni.pdf}
+
+\section{Angular quadrature}
+In double precision, the Lebedev or Beylkin rules suffice.  For higher precision, to exactly integrate angular momenta up to $L$  we employ
+\begin{itemize}
+\item $\theta$: The Gauss-Legendre rule of order $k_\theta = L/2 + 1$ scaled to $[0,\pi]$ ($k$ points can exactly integrate up to $x^{2k-1}$)
+\item $\phi$: $k_\phi = \max(1,2 L)$ equispaced points in $[0,2 \pi)$ ($\phi_i = 2 i \pi / k_\phi$) with weights $2 \pi / k_\phi$.
+ \item The total number of points $k_\theta k_\phi = (l_{\mbox{max}} + 1) \max(4 l_{\mbox{max}},1) $.
+\end{itemize}
+We choose $L = \max(1, 2 l_{\mbox{max}})$, where $l_{\mbox{max}}$ is the maximum angular momentum in the basis, so that we can exactly integrate products.
+
+\begin{table}
+  \begin{center}
+  \caption{Number of angular quadrature points employed for select angular momenta}
+  \begin{tabular}{ccc}
+    $l_{\mbox{max}}$ & Lebedev & GL  \\ \hline
+    0 & 6 & 1 \\
+    1 & 6 & 8 \\
+    2 & 14 & 24 \\
+    3 & 26 & 48 \\
+    4 & 38 & 80 \\
+    5 & 50 & 120 \\
+    6 & 74 & 168 \\ \hline
+  \end{tabular}
+  \end{center}
+\end{table}
+
+\section{Radial quadrature in the b-spline basis}
+
+Between knots, each basis function is just a polynomial of degree $p$.  Since Gauss-Legendre (GL) quadrature with $n$ points can exactly integrate polynomials up to degree $2n-1$, these are the required number of points for various scenarios
+\begin{itemize}
+\item $\int r^l b_i(r) dr$ --- $n \ge (l+p-1)/2$
+\item $\int r^l b_i(r) b_j(r) dr$ --- $n \ge (l+2p-1)/2$
+\end{itemize}
+Pre-tabulated are quadrature points and weights with $n=p+2$, which is adequate for integrals of up to the form $\int r^3 b_i(r) b_j(r)$ which is enough to construct all of the matrices needed to solve atoms.
+
+\section{Projection}
+
+Projecting $f(\vb{r})$ into the representation is accomplished by first integrating over the angular variables with
+
+\begin{eqnarray}
+  \int_0^{2 \pi} d\phi \int_0^\pi d\theta \sin \theta  \N_{lm} (\hat{r}) f(\vb{r}) & = & R_{lm}(r) 
+\end{eqnarray}
+With quadrature points on the unit sphere $\hat{\vb{r}}_\mu$ with
+weights $\omega_\mu$ and radial points $r_\nu$ , this becomes
+\begin{eqnarray}
+  R_{lm}(r_\nu) & = & \sum_\mu \omega_\mu \N_{lm} (\hat{\vb{r}}_\mu) f(r_\nu \hat{\vb{r}}_\mu)
+\end{eqnarray}
+
+The inversion from the function sampled at the grid points can be performed in several ways, paying attention to the ill-conditioning near the origin for the non-zero angular momenta.
+
+By omiting the first (or last) $n$ basis functions the first $n$ derivatives can be forced to be zero on the left (or right).  
+\begin{itemize}
+\item $n=0$ no boundary conditions,
+\item $n=1$ function value is forced to zero,
+\item $n=2$ function value and first derivative are forced to zero,
+\item etc.
+\end{itemize}
+
+\subsection{Weighted least squares and normal equations}
+Adds a penalty to make the problem overall better conditioned and uses the normal equations.  For each $(l,m)$ the LSQ problem is to minimize
+\begin{eqnarray}
+  \sum_\nu  w(r_\nu) \left( R(r_\nu) - \sum_i c_i b_i(r_\nu) \right)^2 + \lambda \sum_i c_i^2
+\end{eqnarray}
+in which the weight is presumably to be chosen as $w(r)=r^2$.
+
+The pentalty term has little effect if $\lambda < \epsilon^2 / \sum_i c_i^2$. 
+
+Setting the variation wrt $c_i$ to zero yields
+\begin{eqnarray}
+  0 & = & - 2 \sum_\nu  w(r_\nu) b_i(r_\nu) \left( R(r_\nu) - \sum_j c_j b_j(r_\nu) \right) + 2\lambda c_i
+\end{eqnarray}
+Tidying up yields
+\begin{eqnarray}
+  \sum_j \left( \sum_\nu w(r_\nu) b_i(r_\nu) b_j(r_\nu) \right) c_j + \lambda c_i & = & \sum_\nu  w(r_\nu) b_i(r_\nu)  R(r_\nu) \\
+  \left( A + \lambda I \right) c & = & b ~ ~\mbox{with} \\
+  a_{ij} & = & \sum_\nu w(r_\nu) b_i(r_\nu) b_j(r_\nu) \\
+  b_i & = & \sum_\nu  w(r_\nu)  b_i(r_\nu) R(r_\nu) 
+\end{eqnarray}
+
+\subsection{Weighted least squares without using normal equations}
+
+We can use a least-squares solver to determine $c$ in
+\begin{eqnarray}
+  w(r_\nu) \sum_i b_i(r_\nu) c_i \approx w(r_\nu) R(r_\nu).
+\end{eqnarray}
+However, it is much more convenient to construct using SVD the pseudo-inverse ($M$) of the matrix
+\begin{eqnarray}
+ a_{\nu i} & = & b_i(r_\nu) w(r_\nu) 
+\end{eqnarray}
+The b-spline basis is strongly linearly independent so there so there should always be a clear divide between zero and non-zero singular values and you can just use the expected number of non-zero singular values.
+
+The fit is obtained directly with a matrix-vector product of $M$ acting on the vector of weighted function values evaluated at the (oversampled) radial points.
+
+\subsection{Projection into the b-spline basis}
+
+Starting from 
+\begin{eqnarray}
+  \sum_j b_j(r) c_j \approx R(r)
+\end{eqnarray}
+we project from the left by $b_i(r)$ to obtain
+\begin{eqnarray}
+  \sum_j \left\langle b_i \middle| b_j\right\rangle c_j & \approx & \left\langle b_i \middle| R \right\rangle
+\end{eqnarray}
+The matrix on the left is just the overlap matrix ($S$) and we expand the integrals on the right using GL quadrature
+\begin{eqnarray}
+\left\langle b_i \middle| R \right\rangle = \sum_\nu b_i(r_\nu) \omega_\nu  R(r_\nu)
+\end{eqnarray}
+Hence, the fitting matrix ($M$) is given by
+\begin{eqnarray}
+   M & = & S^{-1} B
+\end{eqnarray}
+with
+\begin{eqnarray}
+  b_{i \nu} & = & b_i(r_\nu) \omega_\nu 
+\end{eqnarray}
+
+\section{Convolution with an integral operator}
+
+Angular part is diagonal, so focus here on the radial part.
+
+\section{Solid harmonics}
+
+The functions $\N_{lm}(\vb{r})$ are normalized versions of the solid harmonics of Yang ($N_{lm}$).  The unnormalized versions are defined by the recursions (with $m\ge 0$) (the first two are used to recur up the diagonal $m=\pm l$ and the next two are then used to recur $l$ down from $l=m$ to $l=0$ (I inserted commas in the subscripts for clarity)
+\begin{eqnarray}
+  N_{0,0} & = & 1 \\
+  N_{l, l} & = & -\frac{1}{2 l} \left(x N_{l-1,l-1} - y N_{l - 1, -(l - 1)}\right) \\
+  N_{1,-l} & = & -\frac{1}{2 l} \left(y N_{l-1,l-1} + x N_{l - 1, -(l - 1)}\right) \\ 
+  N_{l, m} & = & \frac{1}{(l+m) (l-m)} \left((2 l-1) z N_{l-1, m}-r^2  N_{l-2, m} \right) \\
+  N_{l, -m}& = & \frac{1}{(l+m) (l-m)} \left((2 l-1) z N_{l-1, -m}-r^2 N_{l-2, -m} \right) 
+\end{eqnarray}
+with normalization constants
+\begin{eqnarray}
+  n_{lm} & = & \frac{\sqrt{(2 l +1) \left(l +{| m |}\right)! \left(l -{| m |}\right)!}  2^{-(1 + \delta_{0 m})/2}}{\sqrt{\pi}} .
+\end{eqnarray}
+The normalized harmonics are then
+\begin{eqnarray}
+  \N_{lm}(\vb{r}) & = & n_{lm} N_{lm}(\vb{r}).
+\end{eqnarray}
+The normalization is in the sense that an integral over the unit sphere yields
+\begin{eqnarray}
+  \int_0^{2 \pi} d\phi \int_0^\pi d\theta \sin \theta  \N_{lm} (\hat{r}) \N_{l^\prime m^\prime}(\hat{r})  & = & \delta_{l l^\prime} \delta_{m m^\prime} .
+\end{eqnarray}
+We tabluate the normalization constants and their reciprocals, so we can internally use the simple formulae of the unnormalized harmonics.
+
+Note that $\N_{lm}(\vb{r}) = r^l N_{lm}(\hat{\vb{r}})$.
+
+The first few unnormalized solid harmonics are
+\begin{eqnarray}
+N_{0, 0} & = &  1  \nonumber \\
+N_{1, -1} & = &  -\frac{y}{2}  \nonumber \\
+N_{1, 0} & = &  z  \nonumber \\
+N_{1, 1} & = &  -\frac{x}{2}  \nonumber \\
+N_{2, -2} & = &  \frac{y x}{4}  \nonumber \\
+N_{2, -1} & = &  -\frac{z y}{2}  \nonumber \\
+N_{2, 0} & = &  -\frac{r^{2}}{4}+\frac{3 z^{2}}{4}  \nonumber \\
+N_{2, 1} & = &  -\frac{z x}{2}  \nonumber \\
+N_{2, 2} & = &  \frac{x^{2}}{8}-\frac{y^{2}}{8}  \nonumber \\
+N_{3, -3} & = &  -\frac{1}{16} x^{2} y +\frac{1}{48} y^{3}  \nonumber \\
+N_{3, -2} & = &  \frac{z y x}{4}  \nonumber \\
+N_{3, -1} & = &  \frac{y \left(r^{2}-5 z^{2}\right)}{16}  \nonumber \\
+N_{3, 0} & = &  -\frac{1}{4} r^{2} z +\frac{5}{12} z^{3}  \nonumber \\
+N_{3, 1} & = &  \frac{x \left(r^{2}-5 z^{2}\right)}{16}  \nonumber \\
+N_{3, 2} & = &  \frac{1}{8} z \,x^{2}-\frac{1}{8} z \,y^{2}  \nonumber \\
+N_{3, 3} & = &  -\frac{1}{48} x^{3}+\frac{1}{16} y^{2} x  \nonumber
+\end{eqnarray}
+
+The first few normalized solid harmonics are
+\begin{eqnarray}
+\N_{0, 0} & = &  \frac{1}{2 \sqrt{\pi}}  \nonumber \\
+\N_{1, -1} & = &  -\frac{\sqrt{3} y}{2 \sqrt{\pi}} \nonumber \\
+\N_{1, 0} & = &  \frac{\sqrt{3} z}{2 \sqrt{\pi}} \nonumber \\
+\N_{1, 1} & = &  -\frac{\sqrt{3} x}{2 \sqrt{\pi}} \nonumber \\
+\N_{2, -2} & = &  \frac{\sqrt{15} y x}{2 \sqrt{\pi}} \nonumber \\
+\N_{2, -1} & = &  -\frac{\sqrt{15} z y}{2 \sqrt{\pi}} \nonumber \\
+\N_{2, 0} & = &  -\frac{\sqrt{5} \left(r^{2}-3 z^{2}\right)}{4 \sqrt{\pi}} \nonumber \\
+\N_{2, 1} & = &  -\frac{\sqrt{15} z x}{2 \sqrt{\pi}} \nonumber \\
+\N_{2, 2} & = &  \frac{\sqrt{15} \left(x^{2}-y^{2}\right)}{4 \sqrt{\pi}} \nonumber \\
+\N_{3, -3} & = &  -\frac{3 \sqrt{70} y \left(x^{2}-\frac{y^{2}}{3}\right)}{8 \sqrt{\pi}} \nonumber \\
+\N_{3, -2} & = &  \frac{\sqrt{105} z y x}{2 \sqrt{\pi}} \nonumber \\
+\N_{3, -1} & = &  \frac{y \left(r^{2}-5 z^{2}\right) \sqrt{42}}{8 \sqrt{\pi}} \nonumber \\
+\N_{3, 0} & = &  -\frac{3 \sqrt{7} z \left(r^{2}-\frac{5 z^{2}}{3}\right)}{4 \sqrt{\pi}}  \nonumber \\
+\N_{3, 1} & = &  \frac{x \left(r^{2}-5 z^{2}\right) \sqrt{42}}{8 \sqrt{\pi}} \nonumber \\
+\N_{3, 2} & = &  \frac{\sqrt{105} z \left(x^{2}-y^{2}\right)}{4 \sqrt{\pi}} \nonumber \\
+\N_{3, 3} & = &  -\frac{x \left(x^{2}-3 y^{2}\right) \sqrt{70}}{8 \sqrt{\pi}} \nonumber
+\end{eqnarray}
+
+Derivatives of the unnormalized harmonics are computed as follows with $m\ge 0$, treating as zero $N_{lm}$ with $|m|>l$, and again using commas for clarity
+\begin{eqnarray}
+  \frac{\partial}{\partial x} N_{l,m} & = & \frac{1}{2}\left( N_{l-1,\pm m+1} - N_{l-1,\pm m-1} \right) \nonumber \\
+  \frac{\partial}{\partial y} N_{l,m} & = & \pm \frac{1}{2}\left( N_{l-1,\mp m+1} + N_{l-1,\mp m-1} \right) \nonumber \\
+  \frac{\partial}{\partial z} N_{l,m} & = & N_{l-1,\pm m} . \nonumber
+\end{eqnarray}
+
+\end{document}