|
FLAGLET
1.0b1
Exact wavelets on the ball
|
00001 // FLAGLET package 00002 // Copyright (C) 2012 00003 // Boris Leistedt & Jason McEwen 00004 00005 #include "flaglet.h" 00006 #include <assert.h> 00007 #include <stdlib.h> 00008 #include <stdio.h> 00009 #include <string.h> 00010 #include <math.h> 00011 #include <complex.h> 00012 #include <time.h> 00013 #include <flag.h> 00014 #include <s2let.h> 00015 #include <ssht.h> 00016 00017 void flaglet_random_flmp(complex double *flmp, int L, int P, int seed) 00018 { 00019 int i; 00020 srand( time(NULL) ); 00021 for (i=0; i<L*L*P; i++){ 00022 flmp[i] = (2.0*ran2_dp(seed) - 1.0) + I * (2.0*ran2_dp(seed) - 1.0); 00023 } 00024 } 00025 00026 void flaglet_random_flmp_real(complex double *flmp, int L, int P, int seed) 00027 { 00028 int en, el, m, msign, i, i_op, offset; 00029 int flmsize = ssht_flm_size(L); 00030 00031 for (en=0; en<P; en++) { 00032 offset = en * flmsize; 00033 for (el=0; el<L; el++) { 00034 m = 0; 00035 ssht_sampling_elm2ind(&i, el, m); 00036 flmp[offset+i] = (2.0*ran2_dp(seed) - 1.0); 00037 for (m=1; m<=el; m++) { 00038 ssht_sampling_elm2ind(&i, el, m); 00039 flmp[offset+i] = (2.0*ran2_dp(seed) - 1.0) + I * (2.0*ran2_dp(seed) - 1.0); 00040 ssht_sampling_elm2ind(&i_op, el, -m); 00041 msign = m & 1; 00042 msign = 1 - msign - msign; // (-1)^m 00043 flmp[offset+i_op] = msign * conj(flmp[offset+i]); 00044 } 00045 } 00046 } 00047 } 00048 00049 void flaglet_tiling_test(int B_l, int B_p, int L, int P, int J_min_l, int J_min_p) 00050 { 00051 //double kl, kn; 00052 double *kappa_lp, *kappa0_lp; 00053 flaglet_tiling_axisym_allocate(&kappa_lp, &kappa0_lp, B_l, B_p, L, P); 00054 00055 flaglet_tiling_axisym(kappa_lp, kappa0_lp, B_l, B_p, L, P, J_min_l, J_min_p); 00056 00057 /* 00058 int jl, jp, l, n; 00059 int J_l = s2let_j_max(L, B_l); 00060 int J_p = s2let_j_max(P, B_p); 00061 printf("> KAPPA_0 : "); 00062 for (n = 0; n < P; n++){ 00063 for (l = 0; l < L; l++){ 00064 printf(" %f ,", kappa0_lp[n * L + l]); 00065 } 00066 printf("\n");} 00067 printf("\n"); 00068 for (jp = J_min_p; jp <= J_p; jp++){ 00069 for (jl = J_min_l; jl <= J_l; jl++){ 00070 for (n = 0; n < P; n++){ 00071 for (l = 0; l < L; l++){ 00072 printf(" %f ,",kappa_lp[jp*(J_l+1)*L*P + jl*L*P + n*L + l]); 00073 } 00074 printf("\n");} 00075 printf("\n\n");} 00076 printf("\n\n\n");} 00077 */ 00078 00079 double sum = flaglet_tiling_axisym_check_identity(kappa_lp, kappa0_lp, B_l, B_p, L, P, J_min_l, J_min_p); 00080 printf(" - Identity residuals : %10.5e\n", sum); 00081 00082 free(kappa_lp); 00083 free(kappa0_lp); 00084 } 00085 00086 void flaglet_axisym_wav_lm_test(int B_l, int B_p, int L, int P, int J_min_l, int J_min_p, int seed) 00087 { 00088 clock_t time_start, time_end; 00089 00090 double *wav_lmp, *scal_lmp; 00091 00092 flaglet_axisym_allocate_wav_lmp(&wav_lmp, &scal_lmp, B_l, B_p, L, P); 00093 00094 time_start = clock(); 00095 flaglet_axisym_wav_lmp(wav_lmp, scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p); 00096 time_end = clock(); 00097 printf(" - Generate wavelets : %4.4f seconds\n", 00098 (time_end - time_start) / (double)CLOCKS_PER_SEC); 00099 00100 complex double *f_wav_lmp, *f_scal_lmp, *flmp, *flmp_rec; 00101 flmp = (complex double*)calloc(L * L * P, sizeof(complex double)); 00102 flmp_rec = (complex double*)calloc(L * L * P, sizeof(complex double)); 00103 00104 flaglet_random_flmp(flmp, L, P, seed); 00105 00106 flaglet_axisym_allocate_f_wav_lmp(&f_wav_lmp, &f_scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p); 00107 00108 time_start = clock(); 00109 flaglet_axisym_wav_analysis_lmp(f_wav_lmp, f_scal_lmp, flmp, wav_lmp, scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p); 00110 time_end = clock(); 00111 printf(" - Wavelet analysis : %4.4f seconds\n", 00112 (time_end - time_start) / (double)CLOCKS_PER_SEC); 00113 00114 /* 00115 int n, l, m; 00116 for (n = 0; n < P; n++){ 00117 for (l = 0; l < L; l++){ 00118 for (m = -l; m <= l ; m++){ 00119 int ind = n*L*L + l*l + l + m; 00120 printf("(l,m,n) = (%i,%i,%i) - f_scal_lmp = (%f,%f)\n",l,m,n,creal(f_scal_lmp[ind]),cimag(f_scal_lmp[ind])); 00121 } 00122 } 00123 } 00124 */ 00125 00126 time_start = clock(); 00127 flaglet_axisym_wav_synthesis_lmp(flmp_rec, f_wav_lmp, f_scal_lmp, wav_lmp, scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p); 00128 time_end = clock(); 00129 printf(" - Wavelet synthesis : %4.4f seconds\n", 00130 (time_end - time_start) / (double)CLOCKS_PER_SEC); 00131 00132 /* 00133 for (n = 0; n < P; n++){ 00134 for (l = 0; l < L; l++){ 00135 for (m = -l; m <= l ; m++){ 00136 int ind = n*L*L + l*l + l + m; 00137 if( creal(flmp[ind])-creal(flmp_rec[ind])>0.01 ) 00138 printf("(l,m,n) = (%i,%i,%i) - flmp = (%f,%f) - rec = (%f,%f)\n",l,m,n,creal(flmp[ind]),cimag(flmp[ind]),creal(flmp_rec[ind]),cimag(flmp_rec[ind])); 00139 } 00140 } 00141 } 00142 */ 00143 00144 printf(" - Maximum abs error : %6.5e\n", 00145 maxerr_cplx(flmp, flmp_rec, L * L * P)); 00146 00147 free(flmp); 00148 free(flmp_rec); 00149 free(f_wav_lmp); 00150 free(f_scal_lmp); 00151 free(wav_lmp); 00152 free(scal_lmp); 00153 00154 } 00155 00156 00157 void flaglet_axisym_wav_lm_multires_test(int B_l, int B_p, int L, int P, int J_min_l, int J_min_p, int seed) 00158 { 00159 clock_t time_start, time_end; 00160 00161 double *wav_lmp, *scal_lmp; 00162 00163 flaglet_axisym_allocate_wav_lmp(&wav_lmp, &scal_lmp, B_l, B_p, L, P); 00164 00165 time_start = clock(); 00166 flaglet_axisym_wav_lmp(wav_lmp, scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p); 00167 time_end = clock(); 00168 printf(" - Generate wavelets : %4.4f seconds\n", 00169 (time_end - time_start) / (double)CLOCKS_PER_SEC); 00170 00171 complex double *f_wav_lmp, *f_scal_lmp, *flmp, *flmp_rec; 00172 flmp = (complex double*)calloc(L * L * P, sizeof(complex double)); 00173 flmp_rec = (complex double*)calloc(L * L * P, sizeof(complex double)); 00174 00175 flaglet_random_flmp(flmp, L, P, seed); 00176 00177 flaglet_axisym_allocate_f_wav_multires_lmp(&f_wav_lmp, &f_scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p); 00178 00179 time_start = clock(); 00180 flaglet_axisym_wav_analysis_multires_lmp(f_wav_lmp, f_scal_lmp, flmp, wav_lmp, scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p); 00181 time_end = clock(); 00182 printf(" - Wavelet analysis : %4.4f seconds\n", 00183 (time_end - time_start) / (double)CLOCKS_PER_SEC); 00184 00185 time_start = clock(); 00186 flaglet_axisym_wav_synthesis_multires_lmp(flmp_rec, f_wav_lmp, f_scal_lmp, wav_lmp, scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p); 00187 time_end = clock(); 00188 printf(" - Wavelet synthesis : %4.4f seconds\n", 00189 (time_end - time_start) / (double)CLOCKS_PER_SEC); 00190 00191 printf(" - Maximum abs error : %6.5e\n", 00192 maxerr_cplx(flmp, flmp_rec, L * L * P)); 00193 00194 /* 00195 int n, l, m; 00196 for (n = 0; n < P; n++){ 00197 for (l = 0; l < L; l++){ 00198 for (m = -l; m <= l ; m++){ 00199 int ind = n*L*L + l*l + l + m; 00200 if( creal(flmp[ind])-creal(flmp_rec[ind])>0.01 ) 00201 printf("(l,m,n) = (%i,%i,%i) - flmp = (%f,%f) - rec = (%f,%f)\n",l,m,n,creal(flmp[ind]),cimag(flmp[ind]),creal(flmp_rec[ind]),cimag(flmp_rec[ind])); 00202 } 00203 } 00204 } 00205 */ 00206 00207 free(flmp); 00208 free(flmp_rec); 00209 free(f_wav_lmp); 00210 free(f_scal_lmp); 00211 free(wav_lmp); 00212 free(scal_lmp); 00213 00214 } 00215 00216 00217 void flaglet_axisym_wav_test(double R, int B_l, int B_p, int L, int P, int J_min_l, int J_min_p, int seed) 00218 { 00219 clock_t time_start, time_end; 00220 00221 complex double *f, *f_rec, *flmp, *flmp_rec; 00222 flag_core_allocate_flmn(&flmp, L, P); 00223 flag_core_allocate_flmn(&flmp_rec, L, P); 00224 flag_core_allocate_f(&f, L, P); 00225 flag_core_allocate_f(&f_rec, L, P); 00226 00227 flaglet_random_flmp(flmp, L, P, seed); 00228 00229 double *nodes = (double*)calloc(P, sizeof(double)); 00230 double *weights = (double*)calloc(P, sizeof(double)); 00231 assert(nodes != NULL); 00232 assert(weights != NULL); 00233 flag_spherlaguerre_sampling(nodes, weights, R, P); 00234 00235 flag_core_synthesis(f, flmp, nodes, P, L, P); 00236 00237 free(nodes); 00238 free(weights); 00239 00240 complex double *f_wav, *f_scal; 00241 flaglet_axisym_allocate_f_wav(&f_wav, &f_scal, B_l, B_p, L, P, J_min_l, J_min_p); 00242 00243 time_start = clock(); 00244 flaglet_axisym_wav_analysis(f_wav, f_scal, f, R, B_l, B_p, L, P, J_min_l, J_min_p); 00245 time_end = clock(); 00246 printf(" - Wavelet analysis : %4.4f seconds\n", 00247 (time_end - time_start) / (double)CLOCKS_PER_SEC); 00248 00249 time_start = clock(); 00250 flaglet_axisym_wav_synthesis(f_rec, f_wav, f_scal, R, B_l, B_p, L, P, J_min_l, J_min_p); 00251 time_end = clock(); 00252 printf(" - Wavelet synthesis : %4.4f seconds\n", 00253 (time_end - time_start) / (double)CLOCKS_PER_SEC); 00254 00255 /*printf(" - Maximum abs error : %6.5e\n", 00256 maxerr_cplx(f, f_rec, L*(2*L-1)*P)); 00257 for (n = 0; n < L*(2*L-1)*P; n++){ 00258 //printf("f[%i] = (%f,%f) - rec = (%f,%f)\n",n,creal(f[n]),cimag(f[n]),creal(f_rec[n]),cimag(f_rec[n])); 00259 }*/ 00260 00261 flag_core_analysis(flmp_rec, f_rec, R, L, P); 00262 00263 printf(" - Maximum abs error : %6.5e\n", 00264 maxerr_cplx(flmp, flmp_rec, L*L*P)); 00265 00266 /* 00267 int l, m, n; 00268 for (n = 0; n < P; n++){ 00269 for (l = 0; l < L; l++){ 00270 for (m = -l; m <= l ; m++){ 00271 int ind = n*L*L+l*l+l+m; 00272 printf("(l,m,n) = (%i,%i,%i) - flmp = (%f,%f) - rec = (%f,%f)\n",l,m,n,creal(flmp[ind]),cimag(flmp[ind]),creal(flmp_rec[ind]),cimag(flmp_rec[ind])); 00273 } 00274 } 00275 } 00276 */ 00277 00278 free(f); 00279 free(f_rec); 00280 free(f_wav); 00281 free(f_scal); 00282 } 00283 00284 void flaglet_axisym_wav_real_test(double R, int B_l, int B_p, int L, int P, int J_min_l, int J_min_p, int seed) 00285 { 00286 clock_t time_start, time_end; 00287 00288 complex *flmp, *flmp_rec; 00289 double *f, *f_rec; 00290 flag_core_allocate_flmn(&flmp, L, P); 00291 flag_core_allocate_flmn(&flmp_rec, L, P); 00292 flag_core_allocate_f_real(&f, L, P); 00293 flag_core_allocate_f_real(&f_rec, L, P); 00294 00295 flaglet_random_flmp_real(flmp, L, P, seed); 00296 00297 double *nodes = (double*)calloc(P, sizeof(double)); 00298 double *weights = (double*)calloc(P, sizeof(double)); 00299 assert(nodes != NULL); 00300 assert(weights != NULL); 00301 flag_spherlaguerre_sampling(nodes, weights, R, P); 00302 00303 flag_core_synthesis_real(f, flmp, nodes, P, L, P); 00304 00305 free(nodes); 00306 free(weights); 00307 00308 double *f_wav, *f_scal; 00309 flaglet_axisym_allocate_f_wav_real(&f_wav, &f_scal, B_l, B_p, L, P, J_min_l, J_min_p); 00310 00311 time_start = clock(); 00312 flaglet_axisym_wav_analysis_real(f_wav, f_scal, f, R, B_l, B_p, L, P, J_min_l, J_min_p); 00313 time_end = clock(); 00314 printf(" - Wavelet analysis : %4.4f seconds\n", 00315 (time_end - time_start) / (double)CLOCKS_PER_SEC); 00316 00317 time_start = clock(); 00318 flaglet_axisym_wav_synthesis_real(f_rec, f_wav, f_scal, R, B_l, B_p, L, P, J_min_l, J_min_p); 00319 time_end = clock(); 00320 printf(" - Wavelet synthesis : %4.4f seconds\n", 00321 (time_end - time_start) / (double)CLOCKS_PER_SEC); 00322 00323 /* 00324 printf(" - Maximum abs error : %6.5e\n", 00325 maxerr_cplx(f, f_rec, L*(2*L-1)*P)); 00326 for (n = 0; n < L*(2*L-1)*P; n++){ 00327 printf("f[%i] = (%f,%f) - rec = (%f,%f)\n",n,creal(f[n]),cimag(f[n]),creal(f_rec[n]),cimag(f_rec[n])); 00328 } 00329 */ 00330 00331 flag_core_analysis_real(flmp_rec, f_rec, R, L, P); 00332 00333 printf(" - Maximum abs error : %6.5e\n", 00334 maxerr_cplx(flmp, flmp_rec, L*L*P)); 00335 00336 /* 00337 int l, m, n; 00338 for (n = 0; n < P; n++){ 00339 for (l = 0; l < L; l++){ 00340 for (m = -l; m <= l ; m++){ 00341 int ind = n*L*L+l*l+l+m; 00342 printf("(l,m,n) = (%i,%i,%i) - flmp = (%f,%f) - rec = (%f,%f)\n",l,m,n,creal(flmp[ind]),cimag(flmp[ind]),creal(flmp_rec[ind]),cimag(flmp_rec[ind])); 00343 } 00344 } 00345 } 00346 */ 00347 00348 free(f); 00349 free(f_rec); 00350 free(f_wav); 00351 free(f_scal); 00352 } 00353 00354 00355 void flaglet_axisym_wav_multires_test(double R, int B_l, int B_p, int L, int P, int J_min_l, int J_min_p, int seed) 00356 { 00357 clock_t time_start, time_end; 00358 00359 complex double *f, *f_rec, *flmp, *flmp_rec; 00360 flag_core_allocate_flmn(&flmp, L, P); 00361 flag_core_allocate_flmn(&flmp_rec, L, P); 00362 flag_core_allocate_f(&f, L, P); 00363 flag_core_allocate_f(&f_rec, L, P); 00364 00365 flaglet_random_flmp(flmp, L, P, seed); 00366 00367 double *nodes = (double*)calloc(P, sizeof(double)); 00368 double *weights = (double*)calloc(P, sizeof(double)); 00369 assert(nodes != NULL); 00370 assert(weights != NULL); 00371 flag_spherlaguerre_sampling(nodes, weights, R, P); 00372 00373 flag_core_synthesis(f, flmp, nodes, P, L, P); 00374 00375 free(nodes); 00376 free(weights); 00377 00378 complex double *f_wav, *f_scal; 00379 flaglet_axisym_allocate_f_wav_multires(&f_wav, &f_scal, B_l, B_p, L, P, J_min_l, J_min_p); 00380 00381 time_start = clock(); 00382 flaglet_axisym_wav_analysis_multires(f_wav, f_scal, f, R, B_l, B_p, L, P, J_min_l, J_min_p); 00383 time_end = clock(); 00384 printf(" - Wavelet analysis : %4.4f seconds\n", 00385 (time_end - time_start) / (double)CLOCKS_PER_SEC); 00386 00387 time_start = clock(); 00388 flaglet_axisym_wav_synthesis_multires(f_rec, f_wav, f_scal, R, B_l, B_p, L, P, J_min_l, J_min_p); 00389 time_end = clock(); 00390 printf(" - Wavelet synthesis : %4.4f seconds\n", 00391 (time_end - time_start) / (double)CLOCKS_PER_SEC); 00392 00393 flag_core_analysis(flmp_rec, f_rec, R, L, P); 00394 00395 printf(" - Maximum abs error : %6.5e\n", 00396 maxerr_cplx(flmp, flmp_rec, L*L*P)); 00397 00398 /* 00399 int l, m, n; 00400 for (n = 0; n < P; n++){ 00401 for (l = 0; l < L; l++){ 00402 for (m = -l; m <= l ; m++){ 00403 int ind = n*L*L+l*l+l+m; 00404 if( creal(flmp[ind])-creal(flmp_rec[ind])>0.01 ) 00405 printf("(l,m,n) = (%i,%i,%i) - flmp = (%f,%f) - rec = (%f,%f)\n",l,m,n,creal(flmp[ind]),cimag(flmp[ind]),creal(flmp_rec[ind]),cimag(flmp_rec[ind])); 00406 00407 } 00408 //printf(" %f ", flmp[n*L*L+l*l+l] - flmp_rec[n*L*L+l*l+l]); 00409 } 00410 printf("\n"); 00411 } 00412 */ 00413 00414 free(f); 00415 free(f_rec); 00416 free(f_wav); 00417 free(f_scal); 00418 } 00419 00420 void flaglet_axisym_wav_multires_real_test(double R, int B_l, int B_p, int L, int P, int J_min_l, int J_min_p, int seed) 00421 { 00422 clock_t time_start, time_end; 00423 //int J_l = s2let_j_max(L, B_l); 00424 //int J_p = s2let_j_max(P, B_p); 00425 00426 complex *flmp, *flmp_rec; 00427 double *f, *f_rec; 00428 flag_core_allocate_flmn(&flmp, L, P); 00429 flag_core_allocate_flmn(&flmp_rec, L, P); 00430 flag_core_allocate_f_real(&f, L, P); 00431 flag_core_allocate_f_real(&f_rec, L, P); 00432 00433 flaglet_random_flmp_real(flmp, L, P, seed); 00434 00435 double *nodes = (double*)calloc(P, sizeof(double)); 00436 double *weights = (double*)calloc(P, sizeof(double)); 00437 assert(nodes != NULL); 00438 assert(weights != NULL); 00439 flag_spherlaguerre_sampling(nodes, weights, R, P); 00440 00441 flag_core_synthesis_real(f, flmp, nodes, P, L, P); 00442 00443 free(nodes); 00444 free(weights); 00445 00446 double *f_wav, *f_scal; 00447 flaglet_axisym_allocate_f_wav_multires_real(&f_wav, &f_scal, B_l, B_p, L, P, J_min_l, J_min_p); 00448 00449 time_start = clock(); 00450 flaglet_axisym_wav_analysis_multires_real(f_wav, f_scal, f, R, B_l, B_p, L, P, J_min_l, J_min_p); 00451 time_end = clock(); 00452 printf(" - Wavelet analysis : %4.4f seconds\n", 00453 (time_end - time_start) / (double)CLOCKS_PER_SEC); 00454 00455 time_start = clock(); 00456 flaglet_axisym_wav_synthesis_multires_real(f_rec, f_wav, f_scal, R, B_l, B_p, L, P, J_min_l, J_min_p); 00457 time_end = clock(); 00458 printf(" - Wavelet synthesis : %4.4f seconds\n", 00459 (time_end - time_start) / (double)CLOCKS_PER_SEC); 00460 00461 /* 00462 printf(" - Maximum abs error : %6.5e\n", 00463 maxerr_cplx(f, f_rec, L*(2*L-1)*P)); 00464 for (n = 0; n < L*(2*L-1)*P; n++){ 00465 printf("f[%i] = (%f,%f) - rec = (%f,%f)\n",n,creal(f[n]),cimag(f[n]),creal(f_rec[n]),cimag(f_rec[n])); 00466 } 00467 */ 00468 00469 flag_core_analysis_real(flmp_rec, f_rec, R, L, P); 00470 00471 printf(" - Maximum abs error : %6.5e\n", 00472 maxerr_cplx(flmp, flmp_rec, L*L*P)); 00473 00474 /* 00475 int l, m, n; 00476 for (n = 0; n < P; n++){ 00477 for (l = 0; l < L; l++){ 00478 for (m = -l; m <= l ; m++){ 00479 int ind = n*L*L+l*l+l+m; 00480 if( creal(flmp[ind])-creal(flmp_rec[ind])>0.001 ) 00481 printf("(l,m,n) = (%i,%i,%i) - flmp = (%f,%f) - rec = (%f,%f)\n",l,m,n,creal(flmp[ind]),cimag(flmp[ind]),creal(flmp_rec[ind]),cimag(flmp_rec[ind])); 00482 } 00483 } 00484 } 00485 */ 00486 00487 free(f); 00488 free(f_rec); 00489 free(f_wav); 00490 free(f_scal); 00491 } 00492 00493 00494 void flaglet_axisym_wav_performance_test(double R, int NREPEAT, int NSCALE, int seed, int B_l, int B_p, int J_min_l, int J_min_p) 00495 { 00496 complex double *f, *f_rec, *flmp, *flmp_rec; 00497 clock_t time_start, time_end; 00498 int sc, repeat; 00499 double tottime_analysis = 0, tottime_synthesis = 0; 00500 double accuracy = 0.0; 00501 00502 int L = 2; 00503 int P = 2; 00504 printf(" > B_l = %i - B_p = %i - J_min_l = %i - J_min_p = %i\n",B_l, B_p,J_min_l,J_min_p); 00505 00506 for (sc=0; sc<NSCALE; sc++) { 00507 00508 L *= 2; 00509 P *= 2; 00510 00511 flag_core_allocate_flmn(&flmp, L, P); 00512 00513 double *nodes = (double*)calloc(P, sizeof(double)); 00514 double *weights = (double*)calloc(P, sizeof(double)); 00515 flag_spherlaguerre_sampling(nodes, weights, R, P); 00516 00517 printf("\n > R = %2.2f - L = %i - P = %i \n", R, L, P); 00518 for (repeat=0; repeat<NREPEAT; repeat++){ 00519 00520 //printf(" -> Iteration : %i on %i\n",repeat+1,NREPEAT); 00521 00522 flaglet_random_flmp(flmp, L, P, seed); 00523 00524 flag_core_allocate_f(&f, L, P); 00525 flag_core_allocate_flmn(&flmp_rec, L, P); 00526 00527 time_start = clock(); 00528 flag_core_synthesis(f, flmp, nodes, P, L, P); 00529 time_end = clock(); 00530 //printf(" - (FLAG synthesis : %4.4f seconds)\n", (time_end - time_start) / (double)CLOCKS_PER_SEC); 00531 00532 complex double *f_wav, *f_scal; 00533 flaglet_axisym_allocate_f_wav(&f_wav, &f_scal, B_l, B_p, L, P, J_min_l, J_min_p); 00534 00535 time_start = clock(); 00536 flaglet_axisym_wav_analysis(f_wav, f_scal, f, R, B_l, B_p, L, P, J_min_l, J_min_p); 00537 time_end = clock(); 00538 tottime_analysis += (time_end - time_start) / (double)CLOCKS_PER_SEC; 00539 //printf(" - Duration for wavelet analysis : %4.4f seconds\n", (time_end - time_start) / (double)CLOCKS_PER_SEC); 00540 00541 flag_core_allocate_f(&f_rec, L, P); 00542 00543 time_start = clock(); 00544 flaglet_axisym_wav_synthesis(f_rec, f_wav, f_scal, R, B_l, B_p, L, P, J_min_l, J_min_p); 00545 time_end = clock(); 00546 tottime_synthesis += (time_end - time_start) / (double)CLOCKS_PER_SEC; 00547 //printf(" - Duration for wavelet synthesis : %4.4f seconds\n", (time_end - time_start) / (double)CLOCKS_PER_SEC); 00548 00549 time_start = clock(); 00550 flag_core_analysis(flmp_rec, f_rec, R, L, P); 00551 time_end = clock(); 00552 //printf(" - (FLAG analysis : %4.4f seconds)\n", (time_end - time_start) / (double)CLOCKS_PER_SEC); 00553 00554 accuracy += maxerr_cplx(flmp, flmp_rec, flag_core_flmn_size(L, P)); 00555 //printf(" - Max error on reconstruction : %6.5e\n\n", maxerr_cplx(flmp, flmp_rec, flag_core_flmn_size(L, P))); 00556 00557 free(f); 00558 free(flmp_rec); 00559 free(f_rec); 00560 free(f_wav); 00561 free(f_scal); 00562 00563 } 00564 00565 tottime_analysis = tottime_analysis / (double)NREPEAT; 00566 tottime_synthesis = tottime_synthesis / (double)NREPEAT; 00567 accuracy = accuracy / (double)NREPEAT; 00568 00569 printf(" - Average duration for wavelet analysis : %4.4f seconds\n", tottime_analysis); 00570 printf(" - Average duration for wavelet synthesis : %4.4f seconds\n", tottime_synthesis); 00571 printf(" - Average max error on reconstruction : %6.5e\n", accuracy); 00572 00573 free(flmp); 00574 free(nodes); 00575 free(weights); 00576 00577 } 00578 00579 } 00580 00581 void flaglet_axisym_wav_multires_performance_test(double R, int NREPEAT, int NSCALE, int seed, int B_l, int B_p, int J_min_l, int J_min_p) 00582 { 00583 complex double *f, *f_rec, *flmp, *flmp_rec; 00584 clock_t time_start, time_end; 00585 int sc, repeat; 00586 double tottime_analysis = 0, tottime_synthesis = 0; 00587 double accuracy = 0.0; 00588 00589 int L = 2; 00590 int P = 2; 00591 printf(" > B_l = %i - B_p = %i - J_min_l = %i - J_min_p = %i\n",B_l, B_p,J_min_l,J_min_p); 00592 00593 for (sc=0; sc<NSCALE; sc++) { 00594 00595 L *= 2; 00596 P *= 2; 00597 00598 flag_core_allocate_flmn(&flmp, L, P); 00599 00600 double *nodes = (double*)calloc(P, sizeof(double)); 00601 double *weights = (double*)calloc(P, sizeof(double)); 00602 flag_spherlaguerre_sampling(nodes, weights, R, P); 00603 00604 printf("\n > R = %2.2f - L = %i - P = %i \n", R, L, P); 00605 for (repeat=0; repeat<NREPEAT; repeat++){ 00606 00607 //printf(" -> Iteration : %i on %i\n",repeat+1,NREPEAT); 00608 00609 flaglet_random_flmp(flmp, L, P, seed); 00610 00611 flag_core_allocate_f(&f, L, P); 00612 flag_core_allocate_flmn(&flmp_rec, L, P); 00613 00614 time_start = clock(); 00615 flag_core_synthesis(f, flmp, nodes, P, L, P); 00616 time_end = clock(); 00617 //printf(" - (FLAG synthesis : %4.4f seconds)\n", (time_end - time_start) / (double)CLOCKS_PER_SEC); 00618 00619 complex double *f_wav, *f_scal; 00620 flaglet_axisym_allocate_f_wav_multires(&f_wav, &f_scal, B_l, B_p, L, P, J_min_l, J_min_p); 00621 00622 time_start = clock(); 00623 flaglet_axisym_wav_analysis_multires(f_wav, f_scal, f, R, B_l, B_p, L, P, J_min_l, J_min_p); 00624 time_end = clock(); 00625 tottime_analysis += (time_end - time_start) / (double)CLOCKS_PER_SEC; 00626 //printf(" - Duration for wavelet analysis : %4.4f seconds\n", (time_end - time_start) / (double)CLOCKS_PER_SEC); 00627 00628 flag_core_allocate_f(&f_rec, L, P); 00629 00630 time_start = clock(); 00631 flaglet_axisym_wav_synthesis_multires(f_rec, f_wav, f_scal, R, B_l, B_p, L, P, J_min_l, J_min_p); 00632 time_end = clock(); 00633 tottime_synthesis += (time_end - time_start) / (double)CLOCKS_PER_SEC; 00634 //printf(" - Duration for wavelet synthesis : %4.4f seconds\n", (time_end - time_start) / (double)CLOCKS_PER_SEC); 00635 00636 time_start = clock(); 00637 flag_core_analysis(flmp_rec, f_rec, R, L, P); 00638 time_end = clock(); 00639 //printf(" - (FLAG analysis : %4.4f seconds)\n", (time_end - time_start) / (double)CLOCKS_PER_SEC); 00640 00641 accuracy += maxerr_cplx(flmp, flmp_rec, flag_core_flmn_size(L, P)); 00642 //printf(" - Max error on reconstruction : %6.5e\n\n", maxerr_cplx(flmp, flmp_rec, flag_core_flmn_size(L, P))); 00643 00644 free(f); 00645 free(flmp_rec); 00646 free(f_rec); 00647 free(f_wav); 00648 free(f_scal); 00649 00650 } 00651 00652 tottime_analysis = tottime_analysis / (double)NREPEAT; 00653 tottime_synthesis = tottime_synthesis / (double)NREPEAT; 00654 accuracy = accuracy / (double)NREPEAT; 00655 00656 printf(" - Average duration for wavelet analysis : %4.4f seconds\n", tottime_analysis); 00657 printf(" - Average duration for wavelet synthesis : %4.4f seconds\n", tottime_synthesis); 00658 printf(" - Average max error on reconstruction : %6.5e\n", accuracy); 00659 00660 free(flmp); 00661 free(nodes); 00662 free(weights); 00663 00664 } 00665 00666 } 00667 00668 int main(int argc, char *argv[]) 00669 { 00670 const double R = 1.0; 00671 const int L = 4; 00672 const int P = 4; 00673 const int B_l = 2; 00674 const int B_p = 2; 00675 const int J_min_l = 0; 00676 const int J_min_p = 0; 00677 const int seed = (int)(10000.0*(double)clock()/(double)CLOCKS_PER_SEC); 00678 00679 printf("==========================================================\n"); 00680 printf("PARAMETERS (seed = %i) \n", seed); 00681 printf(" L = %i P = %i Bl = %i Bn = %i Jminl = %i Jminn = %i \n", L, P, B_l, B_p, J_min_l, J_min_p); 00682 printf("----------------------------------------------------------\n"); 00683 printf("> Testing axisymmetric harmonic tiling...\n"); 00684 flaglet_tiling_test(B_l, B_p, L, P, J_min_l, J_min_p); 00685 printf("==========================================================\n"); 00686 printf("> Testing axisymmetric wavelets in harmonics space...\n"); 00687 flaglet_axisym_wav_lm_test(B_l, B_p, L, P, J_min_l, J_min_p, seed); 00688 printf("----------------------------------------------------------\n"); 00689 printf("> Testing multiresolution algorithm in harmonics space...\n"); 00690 flaglet_axisym_wav_lm_multires_test(B_l, B_p, L, P, J_min_l, J_min_p, seed); 00691 printf("==========================================================\n"); 00692 printf("> Testing axisymmetric wavelets in pixel space...\n"); 00693 flaglet_axisym_wav_test(R, B_l, B_p, L, P, J_min_l, J_min_p, seed); 00694 printf("----------------------------------------------------------\n"); 00695 printf("> Testing multiresolution algorithm...\n"); 00696 flaglet_axisym_wav_multires_test(R, B_l, B_p, L, P, J_min_l, J_min_p, seed); 00697 printf("----------------------------------------------------------\n"); 00698 printf("> Testing real axisymmetric wavelets in pixel space...\n"); 00699 flaglet_axisym_wav_real_test(R, B_l, B_p, L, P, J_min_l, J_min_p, seed); 00700 printf("----------------------------------------------------------\n"); 00701 printf("> Testing multiresolution algorithm...\n"); 00702 flaglet_axisym_wav_multires_real_test(R, B_l, B_p, L, P, J_min_l, J_min_p, seed); 00703 /* 00704 const int NREPEAT = 4; 00705 const int NSCALE = 3; 00706 printf("==========================================================\n"); 00707 printf("> Full resolution wavelet transform : performance tests\n"); 00708 flaglet_axisym_wav_performance_test(R, NREPEAT, NSCALE, seed, B_l, B_p, J_min_l, J_min_p); 00709 fflush(NULL); 00710 printf("----------------------------------------------------------\n"); 00711 printf("> Multiresolution wavelet transform : performance tests\n"); 00712 flaglet_axisym_wav_multires_performance_test(R, NREPEAT, NSCALE, seed, B_l, B_p, J_min_l, J_min_p); 00713 fflush(NULL); 00714 */ 00715 printf("==========================================================\n"); 00716 00717 return 0; 00718 }
1.7.6.1
