NFFT  3.3.1
iterS2.c
00001 /*
00002  * Copyright (c) 2002, 2016 Jens Keiner, Stefan Kunis, Daniel Potts
00003  *
00004  * This program is free software; you can redistribute it and/or modify it under
00005  * the terms of the GNU General Public License as published by the Free Software
00006  * Foundation; either version 2 of the License, or (at your option) any later
00007  * version.
00008  *
00009  * This program is distributed in the hope that it will be useful, but WITHOUT
00010  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00011  * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
00012  * details.
00013  *
00014  * You should have received a copy of the GNU General Public License along with
00015  * this program; if not, write to the Free Software Foundation, Inc., 51
00016  * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
00017  */
00018 
00019 /* iterS2 - Iterative reconstruction on the sphere S2 */
00020 
00026 #include "config.h"
00027 
00028 /* Include standard C headers. */
00029 #include <stdio.h>
00030 #include <stdlib.h>
00031 #include <math.h>
00032 #ifdef HAVE_COMPLEX_H
00033 #include <complex.h>
00034 #endif
00035 
00036 /* Include NFFT 3 utilities headers. */
00037 /* Include NFFT3 library header. */
00038 #include "nfft3.h"
00039 #include "infft.h"
00040 
00041 #include "legendre.h"
00042 
00043 static void voronoi_weights_S2(R *w, R *xi, INT M)
00044 {
00045   R *x;
00046   R *y;
00047   R *z;
00048   int j;
00049   int k;
00050   int el;
00051   int Mlocal = M;
00052   int lnew;
00053   int ier;
00054   int *list;
00055   int *lptr;
00056   int *lend;
00057   int *near;
00058   int *next;
00059   R  *dist;
00060   int *ltri;
00061   int *listc;
00062   int nb;
00063   R *xc;
00064   R *yc;
00065   R *zc;
00066   R *rc;
00067   R *vr;
00068   int lp;
00069   int lpl;
00070   int kv;
00071   R a;
00072 
00073   /* Allocate memory for auxilliary arrays. */
00074   x = (R*)X(malloc)(M * sizeof(R));
00075   y = (R*)X(malloc)(M * sizeof(R));
00076   z = (R*)X(malloc)(M * sizeof(R));
00077 
00078   list = (int*)X(malloc)((6*M-12+1)*sizeof(int));
00079   lptr = (int*)X(malloc)((6*M-12+1)*sizeof(int));
00080   lend = (int*)X(malloc)((M+1)*sizeof(int));
00081   near = (int*)X(malloc)((M+1)*sizeof(int));
00082   next = (int*)X(malloc)((M+1)*sizeof(int));
00083   dist = (R*)X(malloc)((M+1)*sizeof(R));
00084   ltri = (int*)X(malloc)((6*M+1)*sizeof(int));
00085   listc = (int*)X(malloc)((6*M-12+1)*sizeof(int));
00086   xc = (R*)X(malloc)((2*M-4+1)*sizeof(R));
00087   yc = (R*)X(malloc)((2*M-4+1)*sizeof(R));
00088   zc = (R*)X(malloc)((2*M-4+1)*sizeof(R));
00089   rc = (R*)X(malloc)((2*M-4+1)*sizeof(R));
00090   vr = (R*)X(malloc)(3*(2*M-4+1)*sizeof(R));
00091 
00092   /* Convert from spherical Coordinates in [0,1/2]x[-1/2,1/2) to Cartesian
00093    * coordinates. */
00094   for (k = 0; k < M; k++)
00095   {
00096     x[k] = SIN(K2PI*xi[2*k+1])*COS(K2PI*xi[2*k]);
00097     y[k] = SIN(K2PI*xi[2*k+1])*SIN(K2PI*xi[2*k]);
00098     z[k] = COS(K2PI*xi[2*k+1]);
00099   }
00100 
00101   /* Generate Delaunay triangulation. */
00102   trmesh_(&Mlocal, x, y, z, list, lptr, lend, &lnew, near, next, dist, &ier);
00103 
00104   /* Check error flag. */
00105   if (ier == 0)
00106   {
00107     /* Get Voronoi vertices. */
00108     crlist_(&Mlocal, &Mlocal, x, y, z, list, lend, lptr, &lnew, ltri, listc, &nb, xc,
00109       yc, zc, rc, &ier);
00110 
00111     if (ier == 0)
00112     {
00113       /* Calcuate sizes of Voronoi regions. */
00114       for (k = 0; k < M; k++)
00115       {
00116         /* Get last neighbour index. */
00117         lpl = lend[k];
00118         lp = lpl;
00119 
00120         j = 0;
00121         vr[3*j] = x[k];
00122         vr[3*j+1] = y[k];
00123         vr[3*j+2] = z[k];
00124 
00125         do
00126         {
00127           j++;
00128           /* Get next neighbour. */
00129           lp = lptr[lp-1];
00130           kv = listc[lp-1];
00131           vr[3*j] = xc[kv-1];
00132           vr[3*j+1] = yc[kv-1];
00133           vr[3*j+2] = zc[kv-1];
00134           /* fprintf(stderr, "lp = %ld\t", lp); */
00135         } while (lp != lpl);
00136 
00137         a = 0;
00138         for (el = 0; el < j; el++)
00139         {
00140           a += areas_(vr, &vr[3*(el+1)],&vr[3*(((el+1)%j)+1)]);
00141         }
00142 
00143         w[k] = a;
00144       }
00145     }
00146   }
00147 
00148   /* Deallocate memory. */
00149   X(free)(x);
00150   X(free)(y);
00151   X(free)(z);
00152 
00153   X(free)(list);
00154   X(free)(lptr);
00155   X(free)(lend);
00156   X(free)(near);
00157   X(free)(next);
00158   X(free)(dist);
00159   X(free)(ltri);
00160   X(free)(listc);
00161   X(free)(xc);
00162   X(free)(yc);
00163   X(free)(zc);
00164   X(free)(rc);
00165   X(free)(vr);
00166 }
00167 
00169 enum boolean {NO = 0, YES = 1};
00170 
00181 int main (int argc, char **argv)
00182 {
00183   int T;
00184   int N;
00185   int M;
00186   int M2;
00187 
00188   int t;                       /* Index variable for testcases                */
00189   nfsft_plan plan;             /* NFSFT plan                                  */
00190   nfsft_plan plan2;            /* NFSFT plan                                  */
00191   solver_plan_complex iplan;           /* NFSFT plan                                  */
00192   int j;                       /*                                             */
00193   int k;                       /*                                             */
00194   int m;                       /*                                             */
00195   int use_nfsft;               /*                                             */
00196   int use_nfft;                /*                                             */
00197   int use_fpt;                 /*                                             */
00198   int cutoff;                  
00199   double threshold;            
00200   double re;
00201   double im;
00202   double a;
00203   double xs;
00204   double *ys;
00205   double *temp;
00206   double _Complex *temp2;
00207   int qlength;
00208   double *qweights;
00209   fftw_plan fplan;
00210   fpt_set set;
00211   int npt;
00212   int npt_exp;
00213   double *alpha, *beta, *gamma;
00214 
00215   /* Read the number of testcases. */
00216   fscanf(stdin,"testcases=%d\n",&T);
00217   fprintf(stderr,"%d\n",T);
00218 
00219   /* Process each testcase. */
00220   for (t = 0; t < T; t++)
00221   {
00222     /* Check if the fast transform shall be used. */
00223     fscanf(stdin,"nfsft=%d\n",&use_nfsft);
00224     fprintf(stderr,"%d\n",use_nfsft);
00225     if (use_nfsft != NO)
00226     {
00227       /* Check if the NFFT shall be used. */
00228       fscanf(stdin,"nfft=%d\n",&use_nfft);
00229       fprintf(stderr,"%d\n",use_nfsft);
00230       if (use_nfft != NO)
00231       {
00232         /* Read the cut-off parameter. */
00233         fscanf(stdin,"cutoff=%d\n",&cutoff);
00234         fprintf(stderr,"%d\n",cutoff);
00235       }
00236       else
00237       {
00238         /* TODO remove this */
00239         /* Initialize unused variable with dummy value. */
00240         cutoff = 1;
00241       }
00242       /* Check if the fast polynomial transform shall be used. */
00243       fscanf(stdin,"fpt=%d\n",&use_fpt);
00244       fprintf(stderr,"%d\n",use_fpt);
00245       if (use_fpt != NO)
00246       {
00247         /* Read the NFSFT threshold parameter. */
00248         fscanf(stdin,"threshold=%lf\n",&threshold);
00249         fprintf(stderr,"%lf\n",threshold);
00250       }
00251       else
00252       {
00253         /* TODO remove this */
00254         /* Initialize unused variable with dummy value. */
00255         threshold = 1000.0;
00256       }
00257     }
00258     else
00259     {
00260       /* TODO remove this */
00261       /* Set dummy values. */
00262       use_nfft = NO;
00263       use_fpt = NO;
00264       cutoff = 3;
00265       threshold = 1000.0;
00266     }
00267 
00268     /* Read the bandwidth. */
00269     fscanf(stdin,"bandwidth=%d\n",&N);
00270     fprintf(stderr,"%d\n",N);
00271 
00272     /* Do precomputation. */
00273     nfsft_precompute(N,threshold,
00274       ((use_nfsft==NO)?(NFSFT_NO_FAST_ALGORITHM):(0U/*NFSFT_NO_DIRECT_ALGORITHM*/)), 0U);
00275 
00276     /* Read the number of nodes. */
00277     fscanf(stdin,"nodes=%d\n",&M);
00278     fprintf(stderr,"%d\n",M);
00279 
00280     /* */
00281     if ((N+1)*(N+1) > M)
00282     {
00283       X(next_power_of_2_exp_int)(N, &npt, &npt_exp);
00284       fprintf(stderr, "npt = %d, npt_exp = %d\n", npt, npt_exp);
00285       fprintf(stderr,"Optimal interpolation!\n");
00286       ys = (double*) nfft_malloc((N+1)*sizeof(double));
00287       temp = (double*) nfft_malloc((2*N+1)*sizeof(double));
00288       temp2 = (double _Complex*) nfft_malloc((N+1)*sizeof(double _Complex));
00289 
00290       a = 0.0;
00291       for (j = 0; j <= N; j++)
00292       {
00293         xs = 2.0 + (2.0*j)/(N+1);
00294         ys[j] = (2.0-((j == 0)?(1.0):(0.0)))*4.0*nfft_bsplines(4,xs);
00295         //fprintf(stdout,"%3d: g(%le) = %le\n",j,xs,ys[j]);
00296         a += ys[j];
00297       }
00298       //fprintf(stdout,"a = %le\n",a);
00299       for (j = 0; j <= N; j++)
00300       {
00301         ys[j] *= 1.0/a;
00302       }
00303 
00304       qlength = 2*N+1;
00305       qweights = (double*) nfft_malloc(qlength*sizeof(double));
00306 
00307       fplan = fftw_plan_r2r_1d(N+1, qweights, qweights, FFTW_REDFT00, 0U);
00308       for (j = 0; j < N+1; j++)
00309       {
00310         qweights[j] = -2.0/(4*j*j-1);
00311       }
00312       fftw_execute(fplan);
00313       qweights[0] *= 0.5;
00314 
00315       for (j = 0; j < N+1; j++)
00316       {
00317         qweights[j] *= 1.0/(2.0*N+1.0);
00318         qweights[2*N+1-1-j] = qweights[j];
00319       }
00320 
00321       fplan = fftw_plan_r2r_1d(2*N+1, temp, temp, FFTW_REDFT00, 0U);
00322       for (j = 0; j <= N; j++)
00323       {
00324         temp[j] = ((j==0 || j == 2*N)?(1.0):(0.5))*ys[j];
00325       }
00326       for (j = N+1; j < 2*N+1; j++)
00327       {
00328         temp[j] = 0.0;
00329       }
00330       fftw_execute(fplan);
00331 
00332       for (j = 0; j < 2*N+1; j++)
00333       {
00334         temp[j] *= qweights[j];
00335       }
00336 
00337       fftw_execute(fplan);
00338 
00339       for (j = 0; j < 2*N+1; j++)
00340       {
00341         temp[j] *= ((j==0 || j == 2*N)?(1.0):(0.5));
00342         if (j <= N)
00343         {
00344           temp2[j] = temp[j];
00345         }
00346       }
00347 
00348       set = fpt_init(1, npt_exp, 0U);
00349 
00350       alpha = (double*) nfft_malloc((N+2)*sizeof(double));
00351       beta = (double*) nfft_malloc((N+2)*sizeof(double));
00352       gamma = (double*) nfft_malloc((N+2)*sizeof(double));
00353 
00354       alpha_al_row(alpha, N, 0);
00355       beta_al_row(beta, N, 0);
00356       gamma_al_row(gamma, N, 0);
00357 
00358       fpt_precompute(set, 0, alpha, beta, gamma, 0, 1000.0);
00359 
00360       fpt_transposed(set,0, temp2, temp2, N, 0U);
00361 
00362       fpt_finalize(set);
00363 
00364       nfft_free(alpha);
00365       nfft_free(beta);
00366       nfft_free(gamma);
00367 
00368       fftw_destroy_plan(fplan);
00369 
00370       nfft_free(qweights);
00371       nfft_free(ys);
00372       nfft_free(temp);
00373     }
00374 
00375     /* Init transform plans. */
00376     nfsft_init_guru(&plan, N, M,
00377       ((use_nfft!=0)?(0U):(NFSFT_USE_NDFT)) |
00378       ((use_fpt!=0)?(0U):(NFSFT_USE_DPT)) | NFSFT_MALLOC_F | NFSFT_MALLOC_X |
00379       NFSFT_MALLOC_F_HAT | NFSFT_NORMALIZED | NFSFT_ZERO_F_HAT,
00380       PRE_PHI_HUT | PRE_PSI | FFTW_INIT |
00381       FFT_OUT_OF_PLACE,
00382       cutoff);
00383 
00384     if ((N+1)*(N+1) > M)
00385     {
00386       solver_init_advanced_complex(&iplan, (nfft_mv_plan_complex*)(&plan), CGNE | PRECOMPUTE_DAMP);
00387     }
00388     else
00389     {
00390       solver_init_advanced_complex(&iplan, (nfft_mv_plan_complex*)(&plan), CGNR | PRECOMPUTE_WEIGHT | PRECOMPUTE_DAMP);
00391     }
00392 
00393     /* Read the nodes and function values. */
00394     for (j = 0; j < M; j++)
00395     {
00396       fscanf(stdin,"%le %le %le %le\n",&plan.x[2*j+1],&plan.x[2*j],&re,&im);
00397       plan.x[2*j+1] = plan.x[2*j+1]/(2.0*KPI);
00398       plan.x[2*j] = plan.x[2*j]/(2.0*KPI);
00399       if (plan.x[2*j] >= 0.5)
00400       {
00401         plan.x[2*j] = plan.x[2*j] - 1;
00402       }
00403       iplan.y[j] = re + _Complex_I * im;
00404       fprintf(stderr,"%le %le %le %le\n",plan.x[2*j+1],plan.x[2*j],
00405         creal(iplan.y[j]),cimag(iplan.y[j]));
00406     }
00407 
00408     /* Read the number of nodes. */
00409     fscanf(stdin,"nodes_eval=%d\n",&M2);
00410     fprintf(stderr,"%d\n",M2);
00411 
00412     /* Init transform plans. */
00413     nfsft_init_guru(&plan2, N, M2,
00414       ((use_nfft!=0)?(0U):(NFSFT_USE_NDFT)) |
00415       ((use_fpt!=0)?(0U):(NFSFT_USE_DPT)) | NFSFT_MALLOC_F | NFSFT_MALLOC_X |
00416       NFSFT_MALLOC_F_HAT | NFSFT_NORMALIZED | NFSFT_ZERO_F_HAT,
00417       PRE_PHI_HUT | PRE_PSI | FFTW_INIT |
00418       FFT_OUT_OF_PLACE,
00419       cutoff);
00420 
00421     /* Read the nodes and function values. */
00422     for (j = 0; j < M2; j++)
00423     {
00424       fscanf(stdin,"%le %le\n",&plan2.x[2*j+1],&plan2.x[2*j]);
00425       plan2.x[2*j+1] = plan2.x[2*j+1]/(2.0*KPI);
00426       plan2.x[2*j] = plan2.x[2*j]/(2.0*KPI);
00427       if (plan2.x[2*j] >= 0.5)
00428       {
00429         plan2.x[2*j] = plan2.x[2*j] - 1;
00430       }
00431       fprintf(stderr,"%le %le\n",plan2.x[2*j+1],plan2.x[2*j]);
00432     }
00433 
00434     nfsft_precompute_x(&plan);
00435 
00436     nfsft_precompute_x(&plan2);
00437 
00438     /* Frequency weights. */
00439     if ((N+1)*(N+1) > M)
00440     {
00441       /* Compute Voronoi weights. */
00442       //voronoi_weights_S2(iplan.w, plan.x, M);
00443 
00444       /* Print out Voronoi weights. */
00445       /*a = 0.0;
00446       for (j = 0; j < plan.M_total; j++)
00447       {
00448         fprintf(stderr,"%le\n",iplan.w[j]);
00449         a += iplan.w[j];
00450       }
00451       fprintf(stderr,"sum = %le\n",a);*/
00452 
00453       for (j = 0; j < plan.N_total; j++)
00454       {
00455         iplan.w_hat[j] = 0.0;
00456       }
00457 
00458       for (k = 0; k <= N; k++)
00459       {
00460         for (j = -k; j <= k; j++)
00461         {
00462           iplan.w_hat[NFSFT_INDEX(k,j,&plan)] = 1.0/(pow(k+1.0,2.0)); /*temp2[j]*/;
00463         }
00464       }
00465     }
00466     else
00467     {
00468       for (j = 0; j < plan.N_total; j++)
00469       {
00470         iplan.w_hat[j] = 0.0;
00471       }
00472 
00473       for (k = 0; k <= N; k++)
00474       {
00475         for (j = -k; j <= k; j++)
00476         {
00477           iplan.w_hat[NFSFT_INDEX(k,j,&plan)] = 1/(pow(k+1.0,2.5));
00478         }
00479       }
00480 
00481       /* Compute Voronoi weights. */
00482       voronoi_weights_S2(iplan.w, plan.x, M);
00483 
00484       /* Print out Voronoi weights. */
00485       a = 0.0;
00486       for (j = 0; j < plan.M_total; j++)
00487       {
00488         fprintf(stderr,"%le\n",iplan.w[j]);
00489         a += iplan.w[j];
00490       }
00491       fprintf(stderr,"sum = %le\n",a);
00492     }
00493 
00494     fprintf(stderr, "N_total = %d\n", plan.N_total);
00495     fprintf(stderr, "M_total = %d\n", plan.M_total);
00496 
00497     /* init some guess */
00498     for (k = 0; k < plan.N_total; k++)
00499     {
00500       iplan.f_hat_iter[k] = 0.0;
00501     }
00502 
00503     /* inverse trafo */
00504     solver_before_loop_complex(&iplan);
00505 
00506     /*for (k = 0; k < plan.M_total; k++)
00507     {
00508       printf("%le %le\n",creal(iplan.r_iter[k]),cimag(iplan.r_iter[k]));
00509     }*/
00510 
00511     for (m = 0; m < 29; m++)
00512     {
00513       fprintf(stderr,"Residual ||r||=%e,\n",sqrt(iplan.dot_r_iter));
00514       solver_loop_one_step_complex(&iplan);
00515     }
00516 
00517     /*CSWAP(iplan.f_hat_iter, plan.f_hat);
00518     nfsft_trafo(&plan);
00519     CSWAP(iplan.f_hat_iter, plan.f_hat);
00520 
00521     a = 0.0;
00522     b = 0.0;
00523     for (k = 0; k < plan.M_total; k++)
00524     {
00525       printf("%le %le %le\n",cabs(iplan.y[k]),cabs(plan.f[k]),
00526         cabs(iplan.y[k]-plan.f[k]));
00527       a += cabs(iplan.y[k]-plan.f[k])*cabs(iplan.y[k]-plan.f[k]);
00528       b += cabs(iplan.y[k])*cabs(iplan.y[k]);
00529     }
00530 
00531     fprintf(stderr,"relative error in 2-norm: %le\n",a/b);*/
00532 
00533     CSWAP(iplan.f_hat_iter, plan2.f_hat);
00534     nfsft_trafo(&plan2);
00535     CSWAP(iplan.f_hat_iter, plan2.f_hat);
00536     for (k = 0; k < plan2.M_total; k++)
00537     {
00538       fprintf(stdout,"%le\n",cabs(plan2.f[k]));
00539     }
00540 
00541     solver_finalize_complex(&iplan);
00542 
00543     nfsft_finalize(&plan);
00544 
00545     nfsft_finalize(&plan2);
00546 
00547     /* Delete precomputed data. */
00548     nfsft_forget();
00549 
00550     if ((N+1)*(N+1) > M)
00551     {
00552       nfft_free(temp2);
00553     }
00554 
00555   } /* Process each testcase. */
00556 
00557   /* Return exit code for successful run. */
00558   return EXIT_SUCCESS;
00559 }
00560 /* \} */