Feat annulus (#75)

* start annulus rotations

* change t to rho

check if the definitions were causing a bug,

start annulus docs

* redo the definitions

* relieve some memory pressure

* plan_ann2cxf

* annulus transform works!

at least for integration

* remove j

* see if the num_threads is causing exit(139)

* add fftw annulus tests

* work on memory safety in annulus

* check modified_jacobi with only polynomial modification (or trivial rational)

* test no modified jacobi

* reduce print format in annulus example

* check only modified_jacobi plan with NULL

* reduce nu

* don't use null

* is it banded_cholfact

* bypass clenshaw too

* delint with -Wall

* return to cholesky

revert all changes to test_transforms_source

* u has length 4

* isolate clang problem to annulus routines

* is it ft_plan_rotannulus?

* caught the bug!

an plan_rotannulus, P->n = XP->n couldn't take place because XP doesn't exist after ft_convert_banded_to_triangular_banded(XP)

* cosmetic fixes to the example and docs

Makefile

@@ -50,6 +50,7 @@ tests:
 
 examples:
 	$(CC) src/ftutilities.c examples/additiontheorem.c $(CFLAGS) -L$(LIBDIR) -l$(LIB) $(LDFLAGS) $(LDLIBS) -o additiontheorem
+	$(CC) src/ftutilities.c examples/annulus.c $(CFLAGS) -L$(LIBDIR) -l$(LIB) $(LDFLAGS) $(LDLIBS) -o annulus
 	$(CC) src/ftutilities.c examples/calculus.c $(CFLAGS) -L$(LIBDIR) -l$(LIB) $(LDFLAGS) $(LDLIBS) -o calculus
 	$(CC) src/ftutilities.c examples/holomorphic.c $(CFLAGS) -L$(LIBDIR) -l$(LIB) $(LDFLAGS) $(LDLIBS) -o holomorphic
 	$(CC) src/ftutilities.c examples/nonlocaldiffusion.c $(CFLAGS) -L$(LIBDIR) -l$(LIB) $(LDFLAGS) $(LDLIBS) -o nonlocaldiffusion
@@ -74,6 +75,7 @@ runtests:
 
 runexamples:
 	./additiontheorem
+	./annulus
 	./calculus
 	./holomorphic
 	./nonlocaldiffusion
@@ -86,6 +88,7 @@ coverage:
 clean:
 	rm -f lib$(LIB).*
 	rm -f additiontheorem
+	rm -f annulus
 	rm -f calculus
 	rm -f holomorphic
 	rm -f nonlocaldiffusion

docs/docs.doc

@@ -246,6 +246,46 @@ f_{2n-2}(r,\theta) = \sum_{\ell=0}^{n-1}\sum_{m=2-2n}^{2n-2} g_{2\ell}^m \left\{
 and \ref ft_execute_cxf2disk converts them back.
 
 
+\subsection ann2cxf
+
+Consider \f$P_n^{t,(\alpha,\beta,\gamma)}(x)\f$ to be orthogonal polynomials in \f$L^2([-1,1], (1-x)^\alpha(1+x)^\beta(t+x)^\gamma\ud x)\f$, for parameter values \f$\{t>1, \alpha,\beta>-1, \gamma\in\mathbb{R}\}\cup\{t=1, \alpha,\beta+\gamma>-1\}\f$.
+
+Real orthonormal annulus polynomials in \f$L^2(\{(r,\theta) : \rho < r < 1, 0 < \theta < 2\pi\}, r^{2\gamma+1}(r^2-\rho^2)^\alpha(1-r^2)^\beta\ud r\ud\theta)\f$ are:
+\f[
+Z_{\ell,m}^{\rho,(\alpha,\beta,\gamma)}(r,\theta) = \sqrt{2} \left(\frac{2}{1-\rho^2}\right)^{\frac{\abs{m}+\alpha+\beta+\gamma+1}{2}} r^{\abs{m}}\tilde{P}_{\frac{\ell-\abs{m}}{2}}^{\frac{1+\rho^2}{1-\rho^2},(\beta,\alpha,\abs{m}+\gamma)}\left(\frac{2r^2-1-\rho^2}{1-\rho^2}\right) \sqrt{\frac{2-\delta_{m,0}}{2\pi}} \left\{\begin{array}{ccc} \cos(m\theta) & {\rm for} & m \ge 0,\\ \sin(\abs{m}\theta) & {\rm for} & m < 0,\end{array}\right.
+\f]
+where the tilde implies that \f$\tilde{P}_n^{t,(\alpha,\beta,\gamma)}(x)\f$ are orthonormal. In the limit as \f$\rho\to0\f$, the annulus polynomials converge to the Zernike polynomials, \f$Z_{\ell,m}^{0,(\alpha,\beta,\gamma)}(r, \theta) = Z_{\ell,m}^{(\alpha+\gamma,\beta)}(r, \theta)\f$. A degree-\f$(2n-2)\f$ expansion in annulus polynomials is given by:
+\f[
+f_{2n-2}(r,\theta) = \sum_{\ell=0}^{2n-2}\sum_{m=-\ell,2}^{+\ell} f_\ell^m Z_{\ell,m}^{\rho,(\alpha,\beta,\gamma)}(r,\theta),
+\f]
+where the \f$,2\f$ in the inner summation index implies that the inner summation runs from \f$m=-\ell\f$ in steps of \f$2\f$ up to \f$+\ell\f$. If annulus polynomial expansion coefficients are organized into the array:
+\f[
+F = \begin{pmatrix}
+f_0^0 & f_1^{-1} & f_1^1 & f_2^{-2} & f_2^2 & \cdots & f_{2n-2}^{2-2n} & f_{2n-2}^{2n-2}\\
+f_2^0 & f_3^{-1} & f_3^1 & f_4^{-2} & f_4^2 & \cdots & 0 & 0\\
+\vdots & \vdots & \vdots &  \vdots &  \vdots & \ddots & \vdots & \vdots\\
+f_{2n-6}^0 & f_{2n-5}^{-1} & f_{2n-5}^1 & f_{2n-4}^{-2} & f_{2n-4}^2 &  & \vdots & \vdots\\
+f_{2n-4}^0 & f_{2n-3}^{-1} & f_{2n-3}^1 & f_{2n-2}^{-2} & f_{2n-2}^2 & \cdots & 0 & 0\\
+f_{2n-2}^0 & 0 & 0 & 0 & 0 & \cdots & 0 & 0\\
+\end{pmatrix},
+\f]
+then \ref ft_plan_ann2cxf creates the appropriate \ref ft_harmonic_plan, and \ref ft_execute_ann2cxf returns the even Chebyshev--Fourier coefficients:
+\f[
+G = \begin{pmatrix}
+g_0^0 & g_0^{-1} & g_0^1 & g_0^{-2} & g_0^2 & \cdots & g_0^{2-2n} & g_0^{2n-2}\\
+g_2^0 & g_2^{-1} & g_2^1 & g_2^{-2} & g_2^2 & \cdots & g_2^{2-2n} & g_2^{2n-2}\\
+\vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\\
+g_{2n-4}^0 & g_{2n-4}^{-1} & g_{2n-4}^1 & g_{2n-4}^{-2} & g_{2n-4}^2 & \cdots & g_{2n-4}^{2-2n} & g_{2n-4}^{2n-2}\\
+g_{2n-2}^0 & 0 & 0 & g_{2n-2}^{-2} & g_{2n-2}^2 & \cdots & g_{2n-2}^{2-2n} & g_{2n-2}^{2n-2}\\
+\end{pmatrix}.
+\f]
+That is:
+\f[
+f_{2n-2}(r,\theta) = \sum_{\ell=0}^{n-1}\sum_{m=2-2n}^{2n-2} g_{2\ell}^m \left\{\begin{array}{ccc} T_{\ell}\left(\frac{2r^2-1-\rho^2}{1-\rho^2}\right) & {\rm for} & m~{\rm even},\\ \sqrt{\frac{2}{1-\rho^2}}rT_{\ell}\left(\frac{2r^2-1-\rho^2}{1-\rho^2}\right) & {\rm for} & m~{\rm odd},\end{array}\right\} \sqrt{\frac{2-\delta_{m,0}}{2\pi}}\left\{\begin{array}{ccc} \cos(m\theta) & {\rm for} & m \ge 0,\\ \sin(\abs{m}\theta) & {\rm for} & m < 0,\end{array}\right.
+\f]
+and \ref ft_execute_cxf2ann converts them back.
+
+
 \subsection rectdisk2cheb
 
 Real orthonormal Dunkl--Xu polynomials in \f$L^2(\mathbb{D}^2, (1-x^2-y^2)^\beta\ud x\ud y)\f$ are:
@@ -436,6 +476,15 @@ r_n & = \cos\left(\frac{2n+1}{4N}\pi\right) = \sin\left(\frac{2N-2n-1}{4N}\pi\ri
 \theta_m & = \frac{2m}{M}\pi,\quad {\rm for} \quad 0 \le m < M.
 \f}
 
+\subsection ft_fftw_annulus_plan
+
+On the annulus, the \f$N\times M\f$ radial--azimuthal grids are:
+\f{align*}
+r_n & = \sqrt{\cos^2\left(\frac{2n+1}{4N}\pi\right) + \rho^2\sin^2\left(\frac{2n+1}{4N}\pi\right)},\\
+& = \sqrt{\sin^2\left(\frac{2N-2n-1}{4N}\pi\right) + \rho^2\cos^2\left(\frac{2N-2n-1}{4N}\pi\right)},\quad {\rm for} \quad 0 \le n < N,\\
+\theta_m & = \frac{2m}{M}\pi,\quad {\rm for} \quad 0 \le m < M.
+\f}
+
 \subsection ft_fftw_rectdisk_plan
 
 On the rectangularized disk, the \f$N\times M\f$ mapped tensor product grids are:

examples/annulus.c

@@ -0,0 +1,102 @@
+#include <fasttransforms.h>
+#include <ftutilities.h>
+
+double f(double x, double y) {return (pow(x, 3))/(x*x+y*y-0.25);};
+
+/*!
+  \example annulus.c
+  In this example, we explore integration of the function:
+  \f[
+  f(x,y) = \frac{x^3}{x^2+y^2-\frac{1}{4}},
+  \f]
+  over the annulus defined by \f$\{(r,\theta) : \frac{2}{3} < r < 1, 0 < \theta < 2\pi\}\f$. We will calculate the integral:
+  \f[
+  \int_0^{2\pi}\int_{\frac{2}{3}}^1 f(r\cos\theta,r\sin\theta)^2r{\rm\,d}r{\rm\,d}\theta,
+  \f]
+  by analyzing the function in an annulus polynomial series.
+*/
+int main(void) {
+    printf("In this example, we explore square integration of a function over \n");
+    printf("the annulus with parameter "MAGENTA("ρ = 2/3")". We analyze the function:\n");
+    printf("\n");
+    printf("\t"MAGENTA("f(x,y) = x³/(x²-y²-1/4)")",\n");
+    printf("\n");
+    printf("on an "MAGENTA("N×M")" tensor product grid defined by:\n");
+    printf("\n");
+    printf("\t"MAGENTA("rₙ = √{cos²[(n+1/2)π/4N] + ρ²sin²[(n+1/2)π/4N]}")", for "MAGENTA("0 ≤ n < N")",\n");
+    printf("\n");
+    printf("and\n");
+    printf("\n");
+    printf("\t"MAGENTA("θₘ = 2π m/M")", for "MAGENTA("0 ≤ m < M")";\n");
+    printf("\n");
+    printf("we convert the function samples to Chebyshev×Fourier coefficients using\n");
+    printf(CYAN("ft_plan_annulus_analysis")" and "CYAN("ft_execute_annulus_analysis")"; and finally, we transform\n");
+    printf("the Chebyshev×Fourier coefficients to annulus harmonic coefficients using\n");
+    printf(CYAN("ft_plan_ann2cxf")" and "CYAN("ft_execute_cxf2ann")".\n");
+    printf("\n");
+    printf("N.B. for the storage pattern of the printed arrays, please consult the\n");
+    printf("documentation. (Arrays are stored in column-major ordering.)\n");
+
+    char * FMT = "%1.3f";
+
+    int N = 8;
+    int M = 4*N-3;
+
+    printf("\n\n"MAGENTA("N = %i")", and "MAGENTA("M = %i")"\n\n", N, M);
+
+    double rho = 2.0/3.0;
+    double r[N], theta[M], F[4*N*N];
+
+    for (int n = 0; n < N; n++) {
+        double t = (N-n-0.5)*M_PI/(2*N);
+        double ct2 = sin(t);
+        double st2 = cos(t);
+        r[n] = sqrt(ct2*ct2+rho*rho*st2*st2);
+    }
+    for (int m = 0; m < M; m++)
+        theta[m] = 2.0*M_PI*m/M;
+
+    printmat("Radial grid "MAGENTA("r"), FMT, r, N, 1);
+    printf("\n");
+    printmat("Azimuthal grid "MAGENTA("θ"), FMT, theta, 1, M);
+    printf("\n");
+
+    for (int m = 0; m < M; m++)
+        for (int n = 0; n < N; n++)
+            F[n+N*m] = f(r[n]*cos(theta[m]), r[n]*sin(theta[m]));
+
+    printf("On the tensor product grid, our function samples are:\n\n");
+
+    printmat("F", FMT, F, N, M);
+    printf("\n");
+
+    double alpha = 0.0, beta = 0.0, gamma = 0.0;
+    ft_harmonic_plan * P = ft_plan_ann2cxf(N, alpha, beta, gamma, rho);
+    ft_annulus_fftw_plan * PA = ft_plan_annulus_analysis(N, M, rho, FT_FFTW_FLAGS);
+
+    ft_execute_annulus_analysis('N', PA, F, N, M);
+    ft_execute_cxf2ann('N', P, F, N, M);
+
+    printf("Its annulus polynomial coefficients are:\n\n");
+
+    printmat("U", FMT, F, N, M);
+    printf("\n");
+
+    printf("The annulus polynomial coefficients are useful for integration.\n");
+    printf("The integral of "MAGENTA("[f(x,y)]^2")" over the annulus is\n");
+    printf("approximately the square of the 2-norm of the coefficients, \n\t");
+    double val = pow(ft_norm_1arg(F, N*M), 2);
+    printf("%1.16f", val);
+    printf(".\n");
+    printf("This compares favourably to the exact result, \n\t");
+    double tval = 5.0*M_PI/8.0*(1675.0/4536.0+9.0*log(3.0)/32.0-3.0*log(7.0)/32.0);
+    printf("%1.16f", tval);
+    printf(".\n");
+    printf("The relative error in the integral is %4.2e.\n", fabs(val-tval)/fabs(tval));
+    printf("This error can be improved upon by increasing "MAGENTA("N")" and "MAGENTA("M")".\n");
+
+    ft_destroy_harmonic_plan(P);
+    ft_destroy_annulus_fftw_plan(PA);
+
+    return 0;
+}

src/banded_source.c

@@ -1502,17 +1502,13 @@ void X(mpsm)(char TRANS, X(modified_plan) * P, FLT * B, int LDB, int N) {
 }
 
 X(modified_plan) * X(plan_modified)(const int n, X(banded) * (*operator_clenshaw)(const int n, const int nc, const FLT * c, const int incc, const X(cop_params) params), const X(cop_params) params, const int nu, const FLT * u, const int nv, const FLT * v, const int verbose) {
+    X(modified_plan) * P = malloc(sizeof(X(modified_plan)));
     if (nv < 1) {
         // polynomial case
         X(banded) * U = operator_clenshaw(n, nu, u, 1, params);
         X(banded_cholfact)(U);
         X(triangular_banded) * R = X(convert_banded_to_triangular_banded)(U);
-        X(modified_plan) * P = malloc(sizeof(X(modified_plan)));
         P->R = R;
-        P->n = n;
-        P->nu = nu;
-        P->nv = nv;
-        return P;
     }
     else {
         // rational case
@@ -1597,14 +1593,13 @@ X(modified_plan) * X(plan_modified)(const int n, X(banded) * (*operator_clenshaw
         X(destroy_banded)(Lt);
         X(destroy_banded)(ULt);
         X(destroy_banded_ql)(F);
-        X(modified_plan) * P = malloc(sizeof(X(modified_plan)));
-        P->n = n;
-        P->nu = nu;
-        P->nv = nv;
         P->K = K;
         P->R = R;
-        return P;
     }
+    P->n = n;
+    P->nu = nu;
+    P->nv = nv;
+    return P;
 }
 
 void Y(execute_jacobi_similarity)(const X(modified_plan) * P, const int n, const FLT * ap, const FLT * bp, FLT * aq, FLT * bq) {