NFFT Logo 3.2.3
iterS2.c
1 /*
2  * Copyright (c) 2002, 2012 Jens Keiner, Stefan Kunis, Daniel Potts
3  *
4  * This program is free software; you can redistribute it and/or modify it under
5  * the terms of the GNU General Public License as published by the Free Software
6  * Foundation; either version 2 of the License, or (at your option) any later
7  * version.
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12  * details.
13  *
14  * You should have received a copy of the GNU General Public License along with
15  * this program; if not, write to the Free Software Foundation, Inc., 51
16  * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  */
18 
19 /* $Id: iterS2.c 3896 2012-10-10 12:19:26Z tovo $ */
20 
21 /* iterS2 - Iterative reconstruction on the sphere S2 */
22 
28 #include "config.h"
29 
30 /* Include standard C headers. */
31 #include <stdio.h>
32 #include <stdlib.h>
33 #include <math.h>
34 #ifdef HAVE_COMPLEX_H
35 #include <complex.h>
36 #endif
37 
38 /* Include NFFT 3 utilities headers. */
39 #include "nfft3util.h"
40 /* Include NFFT3 library header. */
41 #include "nfft3.h"
42 
43 #include "legendre.h"
44 
46 enum boolean {NO = 0, YES = 1};
47 
58 int main (int argc, char **argv)
59 {
60  int T;
61  int N;
62  int M;
63  int M2;
64 
65  int t; /* Index variable for testcases */
66  nfsft_plan plan; /* NFSFT plan */
67  nfsft_plan plan2; /* NFSFT plan */
68  solver_plan_complex iplan; /* NFSFT plan */
69  int j; /* */
70  int k; /* */
71  int m; /* */
72  int use_nfsft; /* */
73  int use_nfft; /* */
74  int use_fpt; /* */
75  int cutoff;
76  double threshold;
77  double re;
78  double im;
79  double a;
80  double *scratch;
81  double xs;
82  double *ys;
83  double *temp;
84  double _Complex *temp2;
85  int qlength;
86  double *qweights;
87  fftw_plan fplan;
88  fpt_set set;
89  int npt;
90  int npt_exp;
91  double *alpha, *beta, *gamma;
92 
93  /* Read the number of testcases. */
94  fscanf(stdin,"testcases=%d\n",&T);
95  fprintf(stderr,"%d\n",T);
96 
97  /* Process each testcase. */
98  for (t = 0; t < T; t++)
99  {
100  /* Check if the fast transform shall be used. */
101  fscanf(stdin,"nfsft=%d\n",&use_nfsft);
102  fprintf(stderr,"%d\n",use_nfsft);
103  if (use_nfsft != NO)
104  {
105  /* Check if the NFFT shall be used. */
106  fscanf(stdin,"nfft=%d\n",&use_nfft);
107  fprintf(stderr,"%d\n",use_nfsft);
108  if (use_nfft != NO)
109  {
110  /* Read the cut-off parameter. */
111  fscanf(stdin,"cutoff=%d\n",&cutoff);
112  fprintf(stderr,"%d\n",cutoff);
113  }
114  else
115  {
116  /* TODO remove this */
117  /* Initialize unused variable with dummy value. */
118  cutoff = 1;
119  }
120  /* Check if the fast polynomial transform shall be used. */
121  fscanf(stdin,"fpt=%d\n",&use_fpt);
122  fprintf(stderr,"%d\n",use_fpt);
123  if (use_fpt != NO)
124  {
125  /* Read the NFSFT threshold parameter. */
126  fscanf(stdin,"threshold=%lf\n",&threshold);
127  fprintf(stderr,"%lf\n",threshold);
128  }
129  else
130  {
131  /* TODO remove this */
132  /* Initialize unused variable with dummy value. */
133  threshold = 1000.0;
134  }
135  }
136  else
137  {
138  /* TODO remove this */
139  /* Set dummy values. */
140  use_nfft = NO;
141  use_fpt = NO;
142  cutoff = 3;
143  threshold = 1000.0;
144  }
145 
146  /* Read the bandwidth. */
147  fscanf(stdin,"bandwidth=%d\n",&N);
148  fprintf(stderr,"%d\n",N);
149 
150  /* Do precomputation. */
151  nfsft_precompute(N,threshold,
152  ((use_nfsft==NO)?(NFSFT_NO_FAST_ALGORITHM):(0U/*NFSFT_NO_DIRECT_ALGORITHM*/)), 0U);
153 
154  /* Read the number of nodes. */
155  fscanf(stdin,"nodes=%d\n",&M);
156  fprintf(stderr,"%d\n",M);
157 
158  /* */
159  if ((N+1)*(N+1) > M)
160  {
161  X(next_power_of_2_exp)(N, &npt, &npt_exp);
162  fprintf(stderr, "npt = %d, npt_exp = %d\n", npt, npt_exp);
163  fprintf(stderr,"Optimal interpolation!\n");
164  scratch = (double*) nfft_malloc(4*sizeof(double));
165  ys = (double*) nfft_malloc((N+1)*sizeof(double));
166  temp = (double*) nfft_malloc((2*N+1)*sizeof(double));
167  temp2 = (double _Complex*) nfft_malloc((N+1)*sizeof(double _Complex));
168 
169  a = 0.0;
170  for (j = 0; j <= N; j++)
171  {
172  xs = 2.0 + (2.0*j)/(N+1);
173  ys[j] = (2.0-((j == 0)?(1.0):(0.0)))*4.0*nfft_bspline(4,xs,scratch);
174  //fprintf(stdout,"%3d: g(%le) = %le\n",j,xs,ys[j]);
175  a += ys[j];
176  }
177  //fprintf(stdout,"a = %le\n",a);
178  for (j = 0; j <= N; j++)
179  {
180  ys[j] *= 1.0/a;
181  }
182 
183  qlength = 2*N+1;
184  qweights = (double*) nfft_malloc(qlength*sizeof(double));
185 
186  fplan = fftw_plan_r2r_1d(N+1, qweights, qweights, FFTW_REDFT00, 0U);
187  for (j = 0; j < N+1; j++)
188  {
189  qweights[j] = -2.0/(4*j*j-1);
190  }
191  fftw_execute(fplan);
192  qweights[0] *= 0.5;
193 
194  for (j = 0; j < N+1; j++)
195  {
196  qweights[j] *= 1.0/(2.0*N+1.0);
197  qweights[2*N+1-1-j] = qweights[j];
198  }
199 
200  fplan = fftw_plan_r2r_1d(2*N+1, temp, temp, FFTW_REDFT00, 0U);
201  for (j = 0; j <= N; j++)
202  {
203  temp[j] = ((j==0 || j == 2*N)?(1.0):(0.5))*ys[j];
204  }
205  for (j = N+1; j < 2*N+1; j++)
206  {
207  temp[j] = 0.0;
208  }
209  fftw_execute(fplan);
210 
211  for (j = 0; j < 2*N+1; j++)
212  {
213  temp[j] *= qweights[j];
214  }
215 
216  fftw_execute(fplan);
217 
218  for (j = 0; j < 2*N+1; j++)
219  {
220  temp[j] *= ((j==0 || j == 2*N)?(1.0):(0.5));
221  if (j <= N)
222  {
223  temp2[j] = temp[j];
224  }
225  }
226 
227  set = fpt_init(1, npt_exp, 0U);
228 
229  alpha = (double*) nfft_malloc((N+2)*sizeof(double));
230  beta = (double*) nfft_malloc((N+2)*sizeof(double));
231  gamma = (double*) nfft_malloc((N+2)*sizeof(double));
232 
233  alpha_al_row(alpha, N, 0);
234  beta_al_row(beta, N, 0);
235  gamma_al_row(gamma, N, 0);
236 
237  fpt_precompute(set, 0, alpha, beta, gamma, 0, 1000.0);
238 
239  fpt_transposed(set,0, temp2, temp2, N, 0U);
240 
241  fpt_finalize(set);
242 
243  nfft_free(alpha);
244  nfft_free(beta);
245  nfft_free(gamma);
246 
247  fftw_destroy_plan(fplan);
248 
249  nfft_free(scratch);
250  nfft_free(qweights);
251  nfft_free(ys);
252  nfft_free(temp);
253  }
254 
255  /* Init transform plans. */
256  nfsft_init_guru(&plan, N, M,
257  ((use_nfft!=0)?(0U):(NFSFT_USE_NDFT)) |
258  ((use_fpt!=0)?(0U):(NFSFT_USE_DPT)) | NFSFT_MALLOC_F | NFSFT_MALLOC_X |
259  NFSFT_MALLOC_F_HAT | NFSFT_NORMALIZED | NFSFT_ZERO_F_HAT,
260  PRE_PHI_HUT | PRE_PSI | FFTW_INIT |
261  FFT_OUT_OF_PLACE,
262  cutoff);
263 
264  if ((N+1)*(N+1) > M)
265  {
266  solver_init_advanced_complex(&iplan, (nfft_mv_plan_complex*)(&plan), CGNE | PRECOMPUTE_DAMP);
267  }
268  else
269  {
270  solver_init_advanced_complex(&iplan, (nfft_mv_plan_complex*)(&plan), CGNR | PRECOMPUTE_WEIGHT | PRECOMPUTE_DAMP);
271  }
272 
273  /* Read the nodes and function values. */
274  for (j = 0; j < M; j++)
275  {
276  fscanf(stdin,"%le %le %le %le\n",&plan.x[2*j+1],&plan.x[2*j],&re,&im);
277  plan.x[2*j+1] = plan.x[2*j+1]/(2.0*PI);
278  plan.x[2*j] = plan.x[2*j]/(2.0*PI);
279  if (plan.x[2*j] >= 0.5)
280  {
281  plan.x[2*j] = plan.x[2*j] - 1;
282  }
283  iplan.y[j] = re + _Complex_I * im;
284  fprintf(stderr,"%le %le %le %le\n",plan.x[2*j+1],plan.x[2*j],
285  creal(iplan.y[j]),cimag(iplan.y[j]));
286  }
287 
288  /* Read the number of nodes. */
289  fscanf(stdin,"nodes_eval=%d\n",&M2);
290  fprintf(stderr,"%d\n",M2);
291 
292  /* Init transform plans. */
293  nfsft_init_guru(&plan2, N, M2,
294  ((use_nfft!=0)?(0U):(NFSFT_USE_NDFT)) |
295  ((use_fpt!=0)?(0U):(NFSFT_USE_DPT)) | NFSFT_MALLOC_F | NFSFT_MALLOC_X |
296  NFSFT_MALLOC_F_HAT | NFSFT_NORMALIZED | NFSFT_ZERO_F_HAT,
297  PRE_PHI_HUT | PRE_PSI | FFTW_INIT |
298  FFT_OUT_OF_PLACE,
299  cutoff);
300 
301  /* Read the nodes and function values. */
302  for (j = 0; j < M2; j++)
303  {
304  fscanf(stdin,"%le %le\n",&plan2.x[2*j+1],&plan2.x[2*j]);
305  plan2.x[2*j+1] = plan2.x[2*j+1]/(2.0*PI);
306  plan2.x[2*j] = plan2.x[2*j]/(2.0*PI);
307  if (plan2.x[2*j] >= 0.5)
308  {
309  plan2.x[2*j] = plan2.x[2*j] - 1;
310  }
311  fprintf(stderr,"%le %le\n",plan2.x[2*j+1],plan2.x[2*j]);
312  }
313 
314  nfsft_precompute_x(&plan);
315 
316  nfsft_precompute_x(&plan2);
317 
318  /* Frequency weights. */
319  if ((N+1)*(N+1) > M)
320  {
321  /* Compute Voronoi weights. */
322  //nfft_voronoi_weights_S2(iplan.w, plan.x, M);
323 
324  /* Print out Voronoi weights. */
325  /*a = 0.0;
326  for (j = 0; j < plan.M_total; j++)
327  {
328  fprintf(stderr,"%le\n",iplan.w[j]);
329  a += iplan.w[j];
330  }
331  fprintf(stderr,"sum = %le\n",a);*/
332 
333  for (j = 0; j < plan.N_total; j++)
334  {
335  iplan.w_hat[j] = 0.0;
336  }
337 
338  for (k = 0; k <= N; k++)
339  {
340  for (j = -k; j <= k; j++)
341  {
342  iplan.w_hat[NFSFT_INDEX(k,j,&plan)] = 1.0/(pow(k+1.0,2.0)); /*temp2[j]*/;
343  }
344  }
345  }
346  else
347  {
348  for (j = 0; j < plan.N_total; j++)
349  {
350  iplan.w_hat[j] = 0.0;
351  }
352 
353  for (k = 0; k <= N; k++)
354  {
355  for (j = -k; j <= k; j++)
356  {
357  iplan.w_hat[NFSFT_INDEX(k,j,&plan)] = 1/(pow(k+1.0,2.5));
358  }
359  }
360 
361  /* Compute Voronoi weights. */
362  nfft_voronoi_weights_S2(iplan.w, plan.x, M);
363 
364  /* Print out Voronoi weights. */
365  a = 0.0;
366  for (j = 0; j < plan.M_total; j++)
367  {
368  fprintf(stderr,"%le\n",iplan.w[j]);
369  a += iplan.w[j];
370  }
371  fprintf(stderr,"sum = %le\n",a);
372  }
373 
374  fprintf(stderr, "N_total = %d\n", plan.N_total);
375  fprintf(stderr, "M_total = %d\n", plan.M_total);
376 
377  /* init some guess */
378  for (k = 0; k < plan.N_total; k++)
379  {
380  iplan.f_hat_iter[k] = 0.0;
381  }
382 
383  /* inverse trafo */
384  solver_before_loop_complex(&iplan);
385 
386  /*for (k = 0; k < plan.M_total; k++)
387  {
388  printf("%le %le\n",creal(iplan.r_iter[k]),cimag(iplan.r_iter[k]));
389  }*/
390 
391  for (m = 0; m < 29; m++)
392  {
393  fprintf(stderr,"Residual ||r||=%e,\n",sqrt(iplan.dot_r_iter));
394  solver_loop_one_step_complex(&iplan);
395  }
396 
397  /*NFFT_SWAP_complex(iplan.f_hat_iter, plan.f_hat);
398  nfsft_trafo(&plan);
399  NFFT_SWAP_complex(iplan.f_hat_iter, plan.f_hat);
400 
401  a = 0.0;
402  b = 0.0;
403  for (k = 0; k < plan.M_total; k++)
404  {
405  printf("%le %le %le\n",cabs(iplan.y[k]),cabs(plan.f[k]),
406  cabs(iplan.y[k]-plan.f[k]));
407  a += cabs(iplan.y[k]-plan.f[k])*cabs(iplan.y[k]-plan.f[k]);
408  b += cabs(iplan.y[k])*cabs(iplan.y[k]);
409  }
410 
411  fprintf(stderr,"relative error in 2-norm: %le\n",a/b);*/
412 
413  NFFT_SWAP_complex(iplan.f_hat_iter, plan2.f_hat);
414  nfsft_trafo(&plan2);
415  NFFT_SWAP_complex(iplan.f_hat_iter, plan2.f_hat);
416  for (k = 0; k < plan2.M_total; k++)
417  {
418  fprintf(stdout,"%le\n",cabs(plan2.f[k]));
419  }
420 
421  solver_finalize_complex(&iplan);
422 
423  nfsft_finalize(&plan);
424 
425  nfsft_finalize(&plan2);
426 
427  /* Delete precomputed data. */
428  nfsft_forget();
429 
430  if ((N+1)*(N+1) > M)
431  {
432  nfft_free(temp2);
433  }
434 
435  } /* Process each testcase. */
436 
437  /* Return exit code for successful run. */
438  return EXIT_SUCCESS;
439 }
440 /* \} */
Generated on Tue Apr 30 2013 by Doxygen 1.8.1