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