NFFT Logo 3.2.2
iterS2.c
00001 /*
00002  * Copyright (c) 2002, 2012 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 /* $Id: iterS2.c 3896 2012-10-10 12:19:26Z tovo $ */
00020 
00021 /* iterS2 - Iterative reconstruction on the sphere S2 */
00022 
00028 #include "config.h"
00029 
00030 /* Include standard C headers. */
00031 #include <stdio.h>
00032 #include <stdlib.h>
00033 #include <math.h>
00034 #ifdef HAVE_COMPLEX_H
00035 #include <complex.h>
00036 #endif
00037 
00038 /* Include NFFT 3 utilities headers. */
00039 #include "nfft3util.h"
00040 /* Include NFFT3 library header. */
00041 #include "nfft3.h"
00042 
00043 #include "legendre.h"
00044 
00046 enum boolean {NO = 0, YES = 1};
00047 
00058 int main (int argc, char **argv)
00059 {
00060   int T;
00061   int N;
00062   int M;
00063   int M2;
00064 
00065   int t;                       /* Index variable for testcases                */
00066   nfsft_plan plan;             /* NFSFT plan                                  */
00067   nfsft_plan plan2;            /* NFSFT plan                                  */
00068   solver_plan_complex iplan;           /* NFSFT plan                                  */
00069   int j;                       /*                                             */
00070   int k;                       /*                                             */
00071   int m;                       /*                                             */
00072   int use_nfsft;               /*                                             */
00073   int use_nfft;                /*                                             */
00074   int use_fpt;                 /*                                             */
00075   int cutoff;                  
00076   double threshold;            
00077   double re;
00078   double im;
00079   double a;
00080   double *scratch;
00081   double xs;
00082   double *ys;
00083   double *temp;
00084   double _Complex *temp2;
00085   int qlength;
00086   double *qweights;
00087   fftw_plan fplan;
00088   fpt_set set;
00089   int npt;
00090   int npt_exp;
00091   double *alpha, *beta, *gamma;
00092 
00093   /* Read the number of testcases. */
00094   fscanf(stdin,"testcases=%d\n",&T);
00095   fprintf(stderr,"%d\n",T);
00096 
00097   /* Process each testcase. */
00098   for (t = 0; t < T; t++)
00099   {
00100     /* Check if the fast transform shall be used. */
00101     fscanf(stdin,"nfsft=%d\n",&use_nfsft);
00102     fprintf(stderr,"%d\n",use_nfsft);
00103     if (use_nfsft != NO)
00104     {
00105       /* Check if the NFFT shall be used. */
00106       fscanf(stdin,"nfft=%d\n",&use_nfft);
00107       fprintf(stderr,"%d\n",use_nfsft);
00108       if (use_nfft != NO)
00109       {
00110         /* Read the cut-off parameter. */
00111         fscanf(stdin,"cutoff=%d\n",&cutoff);
00112         fprintf(stderr,"%d\n",cutoff);
00113       }
00114       else
00115       {
00116         /* TODO remove this */
00117         /* Initialize unused variable with dummy value. */
00118         cutoff = 1;
00119       }
00120       /* Check if the fast polynomial transform shall be used. */
00121       fscanf(stdin,"fpt=%d\n",&use_fpt);
00122       fprintf(stderr,"%d\n",use_fpt);
00123       if (use_fpt != NO)
00124       {
00125         /* Read the NFSFT threshold parameter. */
00126         fscanf(stdin,"threshold=%lf\n",&threshold);
00127         fprintf(stderr,"%lf\n",threshold);
00128       }
00129       else
00130       {
00131         /* TODO remove this */
00132         /* Initialize unused variable with dummy value. */
00133         threshold = 1000.0;
00134       }
00135     }
00136     else
00137     {
00138       /* TODO remove this */
00139       /* Set dummy values. */
00140       use_nfft = NO;
00141       use_fpt = NO;
00142       cutoff = 3;
00143       threshold = 1000.0;
00144     }
00145 
00146     /* Read the bandwidth. */
00147     fscanf(stdin,"bandwidth=%d\n",&N);
00148     fprintf(stderr,"%d\n",N);
00149 
00150     /* Do precomputation. */
00151     nfsft_precompute(N,threshold,
00152       ((use_nfsft==NO)?(NFSFT_NO_FAST_ALGORITHM):(0U/*NFSFT_NO_DIRECT_ALGORITHM*/)), 0U);
00153 
00154     /* Read the number of nodes. */
00155     fscanf(stdin,"nodes=%d\n",&M);
00156     fprintf(stderr,"%d\n",M);
00157 
00158     /* */
00159     if ((N+1)*(N+1) > M)
00160     {
00161       X(next_power_of_2_exp)(N, &npt, &npt_exp);
00162       fprintf(stderr, "npt = %d, npt_exp = %d\n", npt, npt_exp);
00163       fprintf(stderr,"Optimal interpolation!\n");
00164       scratch = (double*) nfft_malloc(4*sizeof(double));
00165       ys = (double*) nfft_malloc((N+1)*sizeof(double));
00166       temp = (double*) nfft_malloc((2*N+1)*sizeof(double));
00167       temp2 = (double _Complex*) nfft_malloc((N+1)*sizeof(double _Complex));
00168 
00169       a = 0.0;
00170       for (j = 0; j <= N; j++)
00171       {
00172         xs = 2.0 + (2.0*j)/(N+1);
00173         ys[j] = (2.0-((j == 0)?(1.0):(0.0)))*4.0*nfft_bspline(4,xs,scratch);
00174         //fprintf(stdout,"%3d: g(%le) = %le\n",j,xs,ys[j]);
00175         a += ys[j];
00176       }
00177       //fprintf(stdout,"a = %le\n",a);
00178       for (j = 0; j <= N; j++)
00179       {
00180         ys[j] *= 1.0/a;
00181       }
00182 
00183       qlength = 2*N+1;
00184       qweights = (double*) nfft_malloc(qlength*sizeof(double));
00185 
00186       fplan = fftw_plan_r2r_1d(N+1, qweights, qweights, FFTW_REDFT00, 0U);
00187       for (j = 0; j < N+1; j++)
00188       {
00189         qweights[j] = -2.0/(4*j*j-1);
00190       }
00191       fftw_execute(fplan);
00192       qweights[0] *= 0.5;
00193 
00194       for (j = 0; j < N+1; j++)
00195       {
00196         qweights[j] *= 1.0/(2.0*N+1.0);
00197         qweights[2*N+1-1-j] = qweights[j];
00198       }
00199 
00200       fplan = fftw_plan_r2r_1d(2*N+1, temp, temp, FFTW_REDFT00, 0U);
00201       for (j = 0; j <= N; j++)
00202       {
00203         temp[j] = ((j==0 || j == 2*N)?(1.0):(0.5))*ys[j];
00204       }
00205       for (j = N+1; j < 2*N+1; j++)
00206       {
00207         temp[j] = 0.0;
00208       }
00209       fftw_execute(fplan);
00210 
00211       for (j = 0; j < 2*N+1; j++)
00212       {
00213         temp[j] *= qweights[j];
00214       }
00215 
00216       fftw_execute(fplan);
00217 
00218       for (j = 0; j < 2*N+1; j++)
00219       {
00220         temp[j] *= ((j==0 || j == 2*N)?(1.0):(0.5));
00221         if (j <= N)
00222         {
00223           temp2[j] = temp[j];
00224         }
00225       }
00226 
00227       set = fpt_init(1, npt_exp, 0U);
00228 
00229       alpha = (double*) nfft_malloc((N+2)*sizeof(double));
00230       beta = (double*) nfft_malloc((N+2)*sizeof(double));
00231       gamma = (double*) nfft_malloc((N+2)*sizeof(double));
00232 
00233       alpha_al_row(alpha, N, 0);
00234       beta_al_row(beta, N, 0);
00235       gamma_al_row(gamma, N, 0);
00236 
00237       fpt_precompute(set, 0, alpha, beta, gamma, 0, 1000.0);
00238 
00239       fpt_transposed(set,0, temp2, temp2, N, 0U);
00240 
00241       fpt_finalize(set);
00242 
00243       nfft_free(alpha);
00244       nfft_free(beta);
00245       nfft_free(gamma);
00246 
00247       fftw_destroy_plan(fplan);
00248 
00249       nfft_free(scratch);
00250       nfft_free(qweights);
00251       nfft_free(ys);
00252       nfft_free(temp);
00253     }
00254 
00255     /* Init transform plans. */
00256     nfsft_init_guru(&plan, N, M,
00257       ((use_nfft!=0)?(0U):(NFSFT_USE_NDFT)) |
00258       ((use_fpt!=0)?(0U):(NFSFT_USE_DPT)) | NFSFT_MALLOC_F | NFSFT_MALLOC_X |
00259       NFSFT_MALLOC_F_HAT | NFSFT_NORMALIZED | NFSFT_ZERO_F_HAT,
00260       PRE_PHI_HUT | PRE_PSI | FFTW_INIT |
00261       FFT_OUT_OF_PLACE,
00262       cutoff);
00263 
00264     if ((N+1)*(N+1) > M)
00265     {
00266       solver_init_advanced_complex(&iplan, (nfft_mv_plan_complex*)(&plan), CGNE | PRECOMPUTE_DAMP);
00267     }
00268     else
00269     {
00270       solver_init_advanced_complex(&iplan, (nfft_mv_plan_complex*)(&plan), CGNR | PRECOMPUTE_WEIGHT | PRECOMPUTE_DAMP);
00271     }
00272 
00273     /* Read the nodes and function values. */
00274     for (j = 0; j < M; j++)
00275     {
00276       fscanf(stdin,"%le %le %le %le\n",&plan.x[2*j+1],&plan.x[2*j],&re,&im);
00277       plan.x[2*j+1] = plan.x[2*j+1]/(2.0*PI);
00278       plan.x[2*j] = plan.x[2*j]/(2.0*PI);
00279       if (plan.x[2*j] >= 0.5)
00280       {
00281         plan.x[2*j] = plan.x[2*j] - 1;
00282       }
00283       iplan.y[j] = re + _Complex_I * im;
00284       fprintf(stderr,"%le %le %le %le\n",plan.x[2*j+1],plan.x[2*j],
00285         creal(iplan.y[j]),cimag(iplan.y[j]));
00286     }
00287 
00288     /* Read the number of nodes. */
00289     fscanf(stdin,"nodes_eval=%d\n",&M2);
00290     fprintf(stderr,"%d\n",M2);
00291 
00292     /* Init transform plans. */
00293     nfsft_init_guru(&plan2, N, M2,
00294       ((use_nfft!=0)?(0U):(NFSFT_USE_NDFT)) |
00295       ((use_fpt!=0)?(0U):(NFSFT_USE_DPT)) | NFSFT_MALLOC_F | NFSFT_MALLOC_X |
00296       NFSFT_MALLOC_F_HAT | NFSFT_NORMALIZED | NFSFT_ZERO_F_HAT,
00297       PRE_PHI_HUT | PRE_PSI | FFTW_INIT |
00298       FFT_OUT_OF_PLACE,
00299       cutoff);
00300 
00301     /* Read the nodes and function values. */
00302     for (j = 0; j < M2; j++)
00303     {
00304       fscanf(stdin,"%le %le\n",&plan2.x[2*j+1],&plan2.x[2*j]);
00305       plan2.x[2*j+1] = plan2.x[2*j+1]/(2.0*PI);
00306       plan2.x[2*j] = plan2.x[2*j]/(2.0*PI);
00307       if (plan2.x[2*j] >= 0.5)
00308       {
00309         plan2.x[2*j] = plan2.x[2*j] - 1;
00310       }
00311       fprintf(stderr,"%le %le\n",plan2.x[2*j+1],plan2.x[2*j]);
00312     }
00313 
00314     nfsft_precompute_x(&plan);
00315 
00316     nfsft_precompute_x(&plan2);
00317 
00318     /* Frequency weights. */
00319     if ((N+1)*(N+1) > M)
00320     {
00321       /* Compute Voronoi weights. */
00322       //nfft_voronoi_weights_S2(iplan.w, plan.x, M);
00323 
00324       /* Print out Voronoi weights. */
00325       /*a = 0.0;
00326       for (j = 0; j < plan.M_total; j++)
00327       {
00328         fprintf(stderr,"%le\n",iplan.w[j]);
00329         a += iplan.w[j];
00330       }
00331       fprintf(stderr,"sum = %le\n",a);*/
00332 
00333       for (j = 0; j < plan.N_total; j++)
00334       {
00335         iplan.w_hat[j] = 0.0;
00336       }
00337 
00338       for (k = 0; k <= N; k++)
00339       {
00340         for (j = -k; j <= k; j++)
00341         {
00342           iplan.w_hat[NFSFT_INDEX(k,j,&plan)] = 1.0/(pow(k+1.0,2.0)); /*temp2[j]*/;
00343         }
00344       }
00345     }
00346     else
00347     {
00348       for (j = 0; j < plan.N_total; j++)
00349       {
00350         iplan.w_hat[j] = 0.0;
00351       }
00352 
00353       for (k = 0; k <= N; k++)
00354       {
00355         for (j = -k; j <= k; j++)
00356         {
00357           iplan.w_hat[NFSFT_INDEX(k,j,&plan)] = 1/(pow(k+1.0,2.5));
00358         }
00359       }
00360 
00361       /* Compute Voronoi weights. */
00362       nfft_voronoi_weights_S2(iplan.w, plan.x, M);
00363 
00364       /* Print out Voronoi weights. */
00365       a = 0.0;
00366       for (j = 0; j < plan.M_total; j++)
00367       {
00368         fprintf(stderr,"%le\n",iplan.w[j]);
00369         a += iplan.w[j];
00370       }
00371       fprintf(stderr,"sum = %le\n",a);
00372     }
00373 
00374     fprintf(stderr, "N_total = %d\n", plan.N_total);
00375     fprintf(stderr, "M_total = %d\n", plan.M_total);
00376 
00377     /* init some guess */
00378     for (k = 0; k < plan.N_total; k++)
00379     {
00380       iplan.f_hat_iter[k] = 0.0;
00381     }
00382 
00383     /* inverse trafo */
00384     solver_before_loop_complex(&iplan);
00385 
00386     /*for (k = 0; k < plan.M_total; k++)
00387     {
00388       printf("%le %le\n",creal(iplan.r_iter[k]),cimag(iplan.r_iter[k]));
00389     }*/
00390 
00391     for (m = 0; m < 29; m++)
00392     {
00393       fprintf(stderr,"Residual ||r||=%e,\n",sqrt(iplan.dot_r_iter));
00394       solver_loop_one_step_complex(&iplan);
00395     }
00396 
00397     /*NFFT_SWAP_complex(iplan.f_hat_iter, plan.f_hat);
00398     nfsft_trafo(&plan);
00399     NFFT_SWAP_complex(iplan.f_hat_iter, plan.f_hat);
00400 
00401     a = 0.0;
00402     b = 0.0;
00403     for (k = 0; k < plan.M_total; k++)
00404     {
00405       printf("%le %le %le\n",cabs(iplan.y[k]),cabs(plan.f[k]),
00406         cabs(iplan.y[k]-plan.f[k]));
00407       a += cabs(iplan.y[k]-plan.f[k])*cabs(iplan.y[k]-plan.f[k]);
00408       b += cabs(iplan.y[k])*cabs(iplan.y[k]);
00409     }
00410 
00411     fprintf(stderr,"relative error in 2-norm: %le\n",a/b);*/
00412 
00413     NFFT_SWAP_complex(iplan.f_hat_iter, plan2.f_hat);
00414     nfsft_trafo(&plan2);
00415     NFFT_SWAP_complex(iplan.f_hat_iter, plan2.f_hat);
00416     for (k = 0; k < plan2.M_total; k++)
00417     {
00418       fprintf(stdout,"%le\n",cabs(plan2.f[k]));
00419     }
00420 
00421     solver_finalize_complex(&iplan);
00422 
00423     nfsft_finalize(&plan);
00424 
00425     nfsft_finalize(&plan2);
00426 
00427     /* Delete precomputed data. */
00428     nfsft_forget();
00429 
00430     if ((N+1)*(N+1) > M)
00431     {
00432       nfft_free(temp2);
00433     }
00434 
00435   } /* Process each testcase. */
00436 
00437   /* Return exit code for successful run. */
00438   return EXIT_SUCCESS;
00439 }
00440 /* \} */

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409