NFFT  3.3.0
iterS2.c
1 /*
2  * Copyright (c) 2002, 2015 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$ */
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 NFFT3 library header. */
40 #include "nfft3.h"
41 #include "infft.h"
42 
43 #include "legendre.h"
44 
45 static void voronoi_weights_S2(R *w, R *xi, INT M)
46 {
47  R *x;
48  R *y;
49  R *z;
50  int j;
51  int k;
52  int el;
53  int Mlocal = M;
54  int lnew;
55  int ier;
56  int *list;
57  int *lptr;
58  int *lend;
59  int *near;
60  int *next;
61  R *dist;
62  int *ltri;
63  int *listc;
64  int nb;
65  R *xc;
66  R *yc;
67  R *zc;
68  R *rc;
69  R *vr;
70  int lp;
71  int lpl;
72  int kv;
73  R a;
74 
75  /* Allocate memory for auxilliary arrays. */
76  x = (R*)X(malloc)(M * sizeof(R));
77  y = (R*)X(malloc)(M * sizeof(R));
78  z = (R*)X(malloc)(M * sizeof(R));
79 
80  list = (int*)X(malloc)((6*M-12+1)*sizeof(int));
81  lptr = (int*)X(malloc)((6*M-12+1)*sizeof(int));
82  lend = (int*)X(malloc)((M+1)*sizeof(int));
83  near = (int*)X(malloc)((M+1)*sizeof(int));
84  next = (int*)X(malloc)((M+1)*sizeof(int));
85  dist = (R*)X(malloc)((M+1)*sizeof(R));
86  ltri = (int*)X(malloc)((6*M+1)*sizeof(int));
87  listc = (int*)X(malloc)((6*M-12+1)*sizeof(int));
88  xc = (R*)X(malloc)((2*M-4+1)*sizeof(R));
89  yc = (R*)X(malloc)((2*M-4+1)*sizeof(R));
90  zc = (R*)X(malloc)((2*M-4+1)*sizeof(R));
91  rc = (R*)X(malloc)((2*M-4+1)*sizeof(R));
92  vr = (R*)X(malloc)(3*(2*M-4+1)*sizeof(R));
93 
94  /* Convert from spherical Coordinates in [0,1/2]x[-1/2,1/2) to Cartesian
95  * coordinates. */
96  for (k = 0; k < M; k++)
97  {
98  x[k] = SIN(K2PI*xi[2*k+1])*COS(K2PI*xi[2*k]);
99  y[k] = SIN(K2PI*xi[2*k+1])*SIN(K2PI*xi[2*k]);
100  z[k] = COS(K2PI*xi[2*k+1]);
101  }
102 
103  /* Generate Delaunay triangulation. */
104  trmesh_(&Mlocal, x, y, z, list, lptr, lend, &lnew, near, next, dist, &ier);
105 
106  /* Check error flag. */
107  if (ier == 0)
108  {
109  /* Get Voronoi vertices. */
110  crlist_(&Mlocal, &Mlocal, x, y, z, list, lend, lptr, &lnew, ltri, listc, &nb, xc,
111  yc, zc, rc, &ier);
112 
113  if (ier == 0)
114  {
115  /* Calcuate sizes of Voronoi regions. */
116  for (k = 0; k < M; k++)
117  {
118  /* Get last neighbour index. */
119  lpl = lend[k];
120  lp = lpl;
121 
122  j = 0;
123  vr[3*j] = x[k];
124  vr[3*j+1] = y[k];
125  vr[3*j+2] = z[k];
126 
127  do
128  {
129  j++;
130  /* Get next neighbour. */
131  lp = lptr[lp-1];
132  kv = listc[lp-1];
133  vr[3*j] = xc[kv-1];
134  vr[3*j+1] = yc[kv-1];
135  vr[3*j+2] = zc[kv-1];
136  /* fprintf(stderr, "lp = %ld\t", lp); */
137  } while (lp != lpl);
138 
139  a = 0;
140  for (el = 0; el < j; el++)
141  {
142  a += areas_(vr, &vr[3*(el+1)],&vr[3*(((el+1)%j)+1)]);
143  }
144 
145  w[k] = a;
146  }
147  }
148  }
149 
150  /* Deallocate memory. */
151  X(free)(x);
152  X(free)(y);
153  X(free)(z);
154 
155  X(free)(list);
156  X(free)(lptr);
157  X(free)(lend);
158  X(free)(near);
159  X(free)(next);
160  X(free)(dist);
161  X(free)(ltri);
162  X(free)(listc);
163  X(free)(xc);
164  X(free)(yc);
165  X(free)(zc);
166  X(free)(rc);
167  X(free)(vr);
168 }
169 
171 enum boolean {NO = 0, YES = 1};
172 
183 int main (int argc, char **argv)
184 {
185  int T;
186  int N;
187  int M;
188  int M2;
189 
190  int t; /* Index variable for testcases */
191  nfsft_plan plan; /* NFSFT plan */
192  nfsft_plan plan2; /* NFSFT plan */
193  solver_plan_complex iplan; /* NFSFT plan */
194  int j; /* */
195  int k; /* */
196  int m; /* */
197  int use_nfsft; /* */
198  int use_nfft; /* */
199  int use_fpt; /* */
200  int cutoff;
201  double threshold;
202  double re;
203  double im;
204  double a;
205  double *scratch;
206  double xs;
207  double *ys;
208  double *temp;
209  double _Complex *temp2;
210  int qlength;
211  double *qweights;
212  fftw_plan fplan;
213  fpt_set set;
214  int npt;
215  int npt_exp;
216  double *alpha, *beta, *gamma;
217 
218  /* Read the number of testcases. */
219  fscanf(stdin,"testcases=%d\n",&T);
220  fprintf(stderr,"%d\n",T);
221 
222  /* Process each testcase. */
223  for (t = 0; t < T; t++)
224  {
225  /* Check if the fast transform shall be used. */
226  fscanf(stdin,"nfsft=%d\n",&use_nfsft);
227  fprintf(stderr,"%d\n",use_nfsft);
228  if (use_nfsft != NO)
229  {
230  /* Check if the NFFT shall be used. */
231  fscanf(stdin,"nfft=%d\n",&use_nfft);
232  fprintf(stderr,"%d\n",use_nfsft);
233  if (use_nfft != NO)
234  {
235  /* Read the cut-off parameter. */
236  fscanf(stdin,"cutoff=%d\n",&cutoff);
237  fprintf(stderr,"%d\n",cutoff);
238  }
239  else
240  {
241  /* TODO remove this */
242  /* Initialize unused variable with dummy value. */
243  cutoff = 1;
244  }
245  /* Check if the fast polynomial transform shall be used. */
246  fscanf(stdin,"fpt=%d\n",&use_fpt);
247  fprintf(stderr,"%d\n",use_fpt);
248  if (use_fpt != NO)
249  {
250  /* Read the NFSFT threshold parameter. */
251  fscanf(stdin,"threshold=%lf\n",&threshold);
252  fprintf(stderr,"%lf\n",threshold);
253  }
254  else
255  {
256  /* TODO remove this */
257  /* Initialize unused variable with dummy value. */
258  threshold = 1000.0;
259  }
260  }
261  else
262  {
263  /* TODO remove this */
264  /* Set dummy values. */
265  use_nfft = NO;
266  use_fpt = NO;
267  cutoff = 3;
268  threshold = 1000.0;
269  }
270 
271  /* Read the bandwidth. */
272  fscanf(stdin,"bandwidth=%d\n",&N);
273  fprintf(stderr,"%d\n",N);
274 
275  /* Do precomputation. */
276  nfsft_precompute(N,threshold,
277  ((use_nfsft==NO)?(NFSFT_NO_FAST_ALGORITHM):(0U/*NFSFT_NO_DIRECT_ALGORITHM*/)), 0U);
278 
279  /* Read the number of nodes. */
280  fscanf(stdin,"nodes=%d\n",&M);
281  fprintf(stderr,"%d\n",M);
282 
283  /* */
284  if ((N+1)*(N+1) > M)
285  {
286  X(next_power_of_2_exp)(N, &npt, &npt_exp);
287  fprintf(stderr, "npt = %d, npt_exp = %d\n", npt, npt_exp);
288  fprintf(stderr,"Optimal interpolation!\n");
289  scratch = (double*) nfft_malloc(4*sizeof(double));
290  ys = (double*) nfft_malloc((N+1)*sizeof(double));
291  temp = (double*) nfft_malloc((2*N+1)*sizeof(double));
292  temp2 = (double _Complex*) nfft_malloc((N+1)*sizeof(double _Complex));
293 
294  a = 0.0;
295  for (j = 0; j <= N; j++)
296  {
297  xs = 2.0 + (2.0*j)/(N+1);
298  ys[j] = (2.0-((j == 0)?(1.0):(0.0)))*4.0*nfft_bspline(4,xs,scratch);
299  //fprintf(stdout,"%3d: g(%le) = %le\n",j,xs,ys[j]);
300  a += ys[j];
301  }
302  //fprintf(stdout,"a = %le\n",a);
303  for (j = 0; j <= N; j++)
304  {
305  ys[j] *= 1.0/a;
306  }
307 
308  qlength = 2*N+1;
309  qweights = (double*) nfft_malloc(qlength*sizeof(double));
310 
311  fplan = fftw_plan_r2r_1d(N+1, qweights, qweights, FFTW_REDFT00, 0U);
312  for (j = 0; j < N+1; j++)
313  {
314  qweights[j] = -2.0/(4*j*j-1);
315  }
316  fftw_execute(fplan);
317  qweights[0] *= 0.5;
318 
319  for (j = 0; j < N+1; j++)
320  {
321  qweights[j] *= 1.0/(2.0*N+1.0);
322  qweights[2*N+1-1-j] = qweights[j];
323  }
324 
325  fplan = fftw_plan_r2r_1d(2*N+1, temp, temp, FFTW_REDFT00, 0U);
326  for (j = 0; j <= N; j++)
327  {
328  temp[j] = ((j==0 || j == 2*N)?(1.0):(0.5))*ys[j];
329  }
330  for (j = N+1; j < 2*N+1; j++)
331  {
332  temp[j] = 0.0;
333  }
334  fftw_execute(fplan);
335 
336  for (j = 0; j < 2*N+1; j++)
337  {
338  temp[j] *= qweights[j];
339  }
340 
341  fftw_execute(fplan);
342 
343  for (j = 0; j < 2*N+1; j++)
344  {
345  temp[j] *= ((j==0 || j == 2*N)?(1.0):(0.5));
346  if (j <= N)
347  {
348  temp2[j] = temp[j];
349  }
350  }
351 
352  set = fpt_init(1, npt_exp, 0U);
353 
354  alpha = (double*) nfft_malloc((N+2)*sizeof(double));
355  beta = (double*) nfft_malloc((N+2)*sizeof(double));
356  gamma = (double*) nfft_malloc((N+2)*sizeof(double));
357 
358  alpha_al_row(alpha, N, 0);
359  beta_al_row(beta, N, 0);
360  gamma_al_row(gamma, N, 0);
361 
362  fpt_precompute(set, 0, alpha, beta, gamma, 0, 1000.0);
363 
364  fpt_transposed(set,0, temp2, temp2, N, 0U);
365 
366  fpt_finalize(set);
367 
368  nfft_free(alpha);
369  nfft_free(beta);
370  nfft_free(gamma);
371 
372  fftw_destroy_plan(fplan);
373 
374  nfft_free(scratch);
375  nfft_free(qweights);
376  nfft_free(ys);
377  nfft_free(temp);
378  }
379 
380  /* Init transform plans. */
381  nfsft_init_guru(&plan, N, M,
382  ((use_nfft!=0)?(0U):(NFSFT_USE_NDFT)) |
383  ((use_fpt!=0)?(0U):(NFSFT_USE_DPT)) | NFSFT_MALLOC_F | NFSFT_MALLOC_X |
384  NFSFT_MALLOC_F_HAT | NFSFT_NORMALIZED | NFSFT_ZERO_F_HAT,
385  PRE_PHI_HUT | PRE_PSI | FFTW_INIT |
386  FFT_OUT_OF_PLACE,
387  cutoff);
388 
389  if ((N+1)*(N+1) > M)
390  {
391  solver_init_advanced_complex(&iplan, (nfft_mv_plan_complex*)(&plan), CGNE | PRECOMPUTE_DAMP);
392  }
393  else
394  {
395  solver_init_advanced_complex(&iplan, (nfft_mv_plan_complex*)(&plan), CGNR | PRECOMPUTE_WEIGHT | PRECOMPUTE_DAMP);
396  }
397 
398  /* Read the nodes and function values. */
399  for (j = 0; j < M; j++)
400  {
401  fscanf(stdin,"%le %le %le %le\n",&plan.x[2*j+1],&plan.x[2*j],&re,&im);
402  plan.x[2*j+1] = plan.x[2*j+1]/(2.0*KPI);
403  plan.x[2*j] = plan.x[2*j]/(2.0*KPI);
404  if (plan.x[2*j] >= 0.5)
405  {
406  plan.x[2*j] = plan.x[2*j] - 1;
407  }
408  iplan.y[j] = re + _Complex_I * im;
409  fprintf(stderr,"%le %le %le %le\n",plan.x[2*j+1],plan.x[2*j],
410  creal(iplan.y[j]),cimag(iplan.y[j]));
411  }
412 
413  /* Read the number of nodes. */
414  fscanf(stdin,"nodes_eval=%d\n",&M2);
415  fprintf(stderr,"%d\n",M2);
416 
417  /* Init transform plans. */
418  nfsft_init_guru(&plan2, N, M2,
419  ((use_nfft!=0)?(0U):(NFSFT_USE_NDFT)) |
420  ((use_fpt!=0)?(0U):(NFSFT_USE_DPT)) | NFSFT_MALLOC_F | NFSFT_MALLOC_X |
421  NFSFT_MALLOC_F_HAT | NFSFT_NORMALIZED | NFSFT_ZERO_F_HAT,
422  PRE_PHI_HUT | PRE_PSI | FFTW_INIT |
423  FFT_OUT_OF_PLACE,
424  cutoff);
425 
426  /* Read the nodes and function values. */
427  for (j = 0; j < M2; j++)
428  {
429  fscanf(stdin,"%le %le\n",&plan2.x[2*j+1],&plan2.x[2*j]);
430  plan2.x[2*j+1] = plan2.x[2*j+1]/(2.0*KPI);
431  plan2.x[2*j] = plan2.x[2*j]/(2.0*KPI);
432  if (plan2.x[2*j] >= 0.5)
433  {
434  plan2.x[2*j] = plan2.x[2*j] - 1;
435  }
436  fprintf(stderr,"%le %le\n",plan2.x[2*j+1],plan2.x[2*j]);
437  }
438 
439  nfsft_precompute_x(&plan);
440 
441  nfsft_precompute_x(&plan2);
442 
443  /* Frequency weights. */
444  if ((N+1)*(N+1) > M)
445  {
446  /* Compute Voronoi weights. */
447  //voronoi_weights_S2(iplan.w, plan.x, M);
448 
449  /* Print out Voronoi weights. */
450  /*a = 0.0;
451  for (j = 0; j < plan.M_total; j++)
452  {
453  fprintf(stderr,"%le\n",iplan.w[j]);
454  a += iplan.w[j];
455  }
456  fprintf(stderr,"sum = %le\n",a);*/
457 
458  for (j = 0; j < plan.N_total; j++)
459  {
460  iplan.w_hat[j] = 0.0;
461  }
462 
463  for (k = 0; k <= N; k++)
464  {
465  for (j = -k; j <= k; j++)
466  {
467  iplan.w_hat[NFSFT_INDEX(k,j,&plan)] = 1.0/(pow(k+1.0,2.0)); /*temp2[j]*/;
468  }
469  }
470  }
471  else
472  {
473  for (j = 0; j < plan.N_total; j++)
474  {
475  iplan.w_hat[j] = 0.0;
476  }
477 
478  for (k = 0; k <= N; k++)
479  {
480  for (j = -k; j <= k; j++)
481  {
482  iplan.w_hat[NFSFT_INDEX(k,j,&plan)] = 1/(pow(k+1.0,2.5));
483  }
484  }
485 
486  /* Compute Voronoi weights. */
487  voronoi_weights_S2(iplan.w, plan.x, M);
488 
489  /* Print out Voronoi weights. */
490  a = 0.0;
491  for (j = 0; j < plan.M_total; j++)
492  {
493  fprintf(stderr,"%le\n",iplan.w[j]);
494  a += iplan.w[j];
495  }
496  fprintf(stderr,"sum = %le\n",a);
497  }
498 
499  fprintf(stderr, "N_total = %d\n", plan.N_total);
500  fprintf(stderr, "M_total = %d\n", plan.M_total);
501 
502  /* init some guess */
503  for (k = 0; k < plan.N_total; k++)
504  {
505  iplan.f_hat_iter[k] = 0.0;
506  }
507 
508  /* inverse trafo */
509  solver_before_loop_complex(&iplan);
510 
511  /*for (k = 0; k < plan.M_total; k++)
512  {
513  printf("%le %le\n",creal(iplan.r_iter[k]),cimag(iplan.r_iter[k]));
514  }*/
515 
516  for (m = 0; m < 29; m++)
517  {
518  fprintf(stderr,"Residual ||r||=%e,\n",sqrt(iplan.dot_r_iter));
519  solver_loop_one_step_complex(&iplan);
520  }
521 
522  /*CSWAP(iplan.f_hat_iter, plan.f_hat);
523  nfsft_trafo(&plan);
524  CSWAP(iplan.f_hat_iter, plan.f_hat);
525 
526  a = 0.0;
527  b = 0.0;
528  for (k = 0; k < plan.M_total; k++)
529  {
530  printf("%le %le %le\n",cabs(iplan.y[k]),cabs(plan.f[k]),
531  cabs(iplan.y[k]-plan.f[k]));
532  a += cabs(iplan.y[k]-plan.f[k])*cabs(iplan.y[k]-plan.f[k]);
533  b += cabs(iplan.y[k])*cabs(iplan.y[k]);
534  }
535 
536  fprintf(stderr,"relative error in 2-norm: %le\n",a/b);*/
537 
538  CSWAP(iplan.f_hat_iter, plan2.f_hat);
539  nfsft_trafo(&plan2);
540  CSWAP(iplan.f_hat_iter, plan2.f_hat);
541  for (k = 0; k < plan2.M_total; k++)
542  {
543  fprintf(stdout,"%le\n",cabs(plan2.f[k]));
544  }
545 
546  solver_finalize_complex(&iplan);
547 
548  nfsft_finalize(&plan);
549 
550  nfsft_finalize(&plan2);
551 
552  /* Delete precomputed data. */
553  nfsft_forget();
554 
555  if ((N+1)*(N+1) > M)
556  {
557  nfft_free(temp2);
558  }
559 
560  } /* Process each testcase. */
561 
562  /* Return exit code for successful run. */
563  return EXIT_SUCCESS;
564 }
565 /* \} */
double * x
the nodes for ,
Definition: nfft3.h:576
int main(int argc, char **argv)
The main program.
Definition: iterS2.c:183
fftw_complex * f_hat
Fourier coefficients.
Definition: nfft3.h:576
double * w
weighting factors
Definition: nfft3.h:786
data structure for an NFSFT (nonequispaced fast spherical Fourier transform) plan with double precisi...
Definition: nfft3.h:576
double dot_r_iter
weighted dotproduct of r_iter
Definition: nfft3.h:786
Holds data for a set of cascade summations.
Definition: fpt.c:96
void nfft_free(void *p)
NFFT_INT M_total
Total number of samples.
Definition: nfft3.h:576
#define X(name)
Include header for C99 complex datatype.
Definition: fastsum.h:53
void * nfft_malloc(size_t n)
NFFT_INT N_total
Total number of Fourier coefficients.
Definition: nfft3.h:576
fftw_complex * f
Samples.
Definition: nfft3.h:576
fftw_complex * y
right hand side, samples
Definition: nfft3.h:786
data structure for an inverse NFFT plan with double precision
Definition: nfft3.h:786
#define CSWAP(x, y)
Swap two vectors.
Definition: infft.h:152
double * w_hat
damping factors
Definition: nfft3.h:786
fftw_complex * f_hat_iter
iterative solution
Definition: nfft3.h:786