![]() |
NFFT
3.3.1
|
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 /* \} */