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#include <stdio.h> 20#include <math.h> 21#include <string.h> 22#include <stdlib.h> 23#include <complex.h> 24 25#define @NFFT_PRECISION_MACRO@ 26 27#include "nfft3mp.h" 28 29/** 30 * \defgroup examples_solver_glacier Reconstruction of a glacier from \ 31 scattered data 32 * \ingroup examples_solver 33 * \{ 34 */ 35 36/** Generalised Sobolev weight */ 37static NFFT_R my_weight(NFFT_R z, NFFT_R a, NFFT_R b, NFFT_R c) 38{ 39 return NFFT_M(pow)(NFFT_K(0.25) - z * z, b) / (c + NFFT_M(pow)(NFFT_M(fabs)(z), NFFT_K(2.0) * a)); 40} 41 42/** Reconstruction routine */ 43static void glacier(int N, int M) 44{ 45 int j, k, k0, k1, l, my_N[2], my_n[2]; 46 NFFT_R tmp_y; 47 NFFT(plan) p; 48 SOLVER(plan_complex) ip; 49 FILE* fp; 50 51 /* initialise p */ 52 my_N[0] = N; 53 my_n[0] = (int)NFFT(next_power_of_2)(N); 54 my_N[1] = N; 55 my_n[1] = (int)NFFT(next_power_of_2)(N); 56 NFFT(init_guru)(&p, 2, my_N, M, my_n, 6, 57 PRE_PHI_HUT | PRE_FULL_PSI | 58 MALLOC_X | MALLOC_F_HAT | MALLOC_F | 59 FFTW_INIT | FFT_OUT_OF_PLACE, 60 FFTW_MEASURE | FFTW_DESTROY_INPUT); 61 62 /* initialise ip, specific */ 63 SOLVER(init_advanced_complex)(&ip, (NFFT(mv_plan_complex)*) (&p), 64 CGNE | PRECOMPUTE_DAMP); 65 fprintf(stderr, "Using the generic solver!"); 66 67 /* init nodes */ 68 fp = fopen("input_data.dat", "r"); 69 for (j = 0; j < p.M_total; j++) 70 { 71 fscanf(fp, "%" NFFT__FES__ " %" NFFT__FES__ " %" NFFT__FES__, &p.x[2 * j + 0], &p.x[2 * j + 1], &tmp_y); 72 ip.y[j] = tmp_y; 73 } 74 fclose(fp); 75 76 /* precompute psi */ 77 if (p.flags & PRE_ONE_PSI) 78 NFFT(precompute_one_psi)(&p); 79 80 /* initialise damping factors */ 81 if (ip.flags & PRECOMPUTE_DAMP) 82 for (k0 = 0; k0 < p.N[0]; k0++) 83 for (k1 = 0; k1 < p.N[1]; k1++) 84 ip.w_hat[k0 * p.N[1] + k1] = my_weight(((NFFT_R) ((NFFT_R)(k0) - (NFFT_R)(p.N[0]) / NFFT_K(2.0))) / ((NFFT_R)(p.N[0])), 85 NFFT_K(0.5), NFFT_K(3.0), NFFT_K(0.001)) 86 * my_weight((((NFFT_R)(k1) - (NFFT_R)(p.N[1]) / NFFT_K(2.0))) / ((NFFT_R)(p.N[1])), NFFT_K(0.5), NFFT_K(3.0), NFFT_K(0.001)); 87 88 /* init some guess */ 89 for (k = 0; k < p.N_total; k++) 90 ip.f_hat_iter[k] = NFFT_K(0.0); 91 92 /* inverse trafo */ 93 SOLVER(before_loop_complex)(&ip); 94 95 for (l = 0; l < 40; l++) 96 { 97 fprintf(stderr, "Residual ||r||=%" NFFT__FES__ ",\n", NFFT_M(sqrt)(ip.dot_r_iter)); 98 SOLVER(loop_one_step_complex)(&ip); 99 } 100 101 for (k = 0; k < p.N_total; k++) 102 printf("%" NFFT__FES__ " %" NFFT__FES__ "\n", NFFT_M(creal)(ip.f_hat_iter[k]), NFFT_M(cimag)(ip.f_hat_iter[k])); 103 104 SOLVER(finalize_complex)(&ip); 105 NFFT(finalize)(&p); 106} 107 108/** Reconstruction routine with cross validation */ 109static void glacier_cv(int N, int M, int M_cv, unsigned solver_flags) 110{ 111 int j, k, k0, k1, l, my_N[2], my_n[2]; 112 NFFT_R tmp_y, r; 113 NFFT(plan) p, cp; 114 SOLVER(plan_complex) ip; 115 NFFT_C* cp_y; 116 FILE* fp; 117 int M_re = M - M_cv; 118 119 /* initialise p for reconstruction */ 120 my_N[0] = N; 121 my_n[0] = (int)NFFT(next_power_of_2)(N); 122 my_N[1] = N; 123 my_n[1] = (int)NFFT(next_power_of_2)(N); 124 NFFT(init_guru)(&p, 2, my_N, M_re, my_n, 6, 125 PRE_PHI_HUT | PRE_FULL_PSI | 126 MALLOC_X | MALLOC_F_HAT | MALLOC_F | 127 FFTW_INIT | FFT_OUT_OF_PLACE, 128 FFTW_MEASURE | FFTW_DESTROY_INPUT); 129 130 /* initialise ip, specific */ 131 SOLVER(init_advanced_complex)(&ip, (NFFT(mv_plan_complex)*) (&p), solver_flags); 132 133 /* initialise cp for validation */ 134 cp_y = (NFFT_C*) NFFT(malloc)((size_t)(M) * sizeof(NFFT_C)); 135 NFFT(init_guru)(&cp, 2, my_N, M, my_n, 6, 136 PRE_PHI_HUT | PRE_FULL_PSI | 137 MALLOC_X | MALLOC_F | 138 FFTW_INIT | FFT_OUT_OF_PLACE, 139 FFTW_MEASURE | FFTW_DESTROY_INPUT); 140 141 cp.f_hat = ip.f_hat_iter; 142 143 /* set up data in cp and cp_y */ 144 fp = fopen("input_data.dat", "r"); 145 for (j = 0; j < cp.M_total; j++) 146 { 147 fscanf(fp, "%" NFFT__FES__ " %" NFFT__FES__ " %" NFFT__FES__, &cp.x[2 * j + 0], &cp.x[2 * j + 1], &tmp_y); 148 cp_y[j] = tmp_y; 149 } 150 fclose(fp); 151 152 /* copy part of the data to p and ip */ 153 for (j = 0; j < p.M_total; j++) 154 { 155 p.x[2 * j + 0] = cp.x[2 * j + 0]; 156 p.x[2 * j + 1] = cp.x[2 * j + 1]; 157 ip.y[j] = tmp_y; 158 } 159 160 /* precompute psi */ 161 if (p.flags & PRE_ONE_PSI) 162 NFFT(precompute_one_psi)(&p); 163 164 /* precompute psi */ 165 if (cp.flags & PRE_ONE_PSI) 166 NFFT(precompute_one_psi)(&cp); 167 168 /* initialise damping factors */ 169 if (ip.flags & PRECOMPUTE_DAMP) 170 for (k0 = 0; k0 < p.N[0]; k0++) 171 for (k1 = 0; k1 < p.N[1]; k1++) 172 ip.w_hat[k0 * p.N[1] + k1] = my_weight((((NFFT_R)(k0) - (NFFT_R)(p.N[0]) / NFFT_K(2.0))) / ((NFFT_R)(p.N[0])), 173 NFFT_K(0.5), NFFT_K(3.0), NFFT_K(0.001)) 174 * my_weight((((NFFT_R)(k1) - (NFFT_R)(p.N[1]) / NFFT_K(2.0))) / ((NFFT_R)(p.N[1])), NFFT_K(0.5), NFFT_K(3.0), NFFT_K(0.001)); 175 176 /* init some guess */ 177 for (k = 0; k < p.N_total; k++) 178 ip.f_hat_iter[k] = NFFT_K(0.0); 179 180 /* inverse trafo */ 181 SOLVER(before_loop_complex)(&ip); 182 // fprintf(stderr,"iteration starts,\t"); 183 for (l = 0; l < 40; l++) 184 SOLVER(loop_one_step_complex)(&ip); 185 186 //fprintf(stderr,"r=%1.2e, ",sqrt(ip.dot_r_iter)/M_re); 187 188 NFFT_CSWAP(p.f_hat, ip.f_hat_iter); 189 NFFT(trafo)(&p); 190 NFFT_CSWAP(p.f_hat, ip.f_hat_iter); 191 NFFT(upd_axpy_complex)(p.f, -1, ip.y, M_re); 192 r = NFFT_M(sqrt)(NFFT(dot_complex)(p.f, M_re) / NFFT(dot_complex)(cp_y, M)); 193 fprintf(stderr, "r=%1.2" NFFT__FES__ ", ", r); 194 printf("$%1.1" NFFT__FES__ "$ & ", r); 195 196 NFFT(trafo)(&cp); 197 NFFT(upd_axpy_complex)(&cp.f[M_re], -1, &cp_y[M_re], M_cv); 198 r = NFFT_M(sqrt)(NFFT(dot_complex)(&cp.f[M_re], M_cv) / NFFT(dot_complex)(cp_y, M)); 199 fprintf(stderr, "r_1=%1.2" NFFT__FES__ "\t", r); 200 printf("$%1.1" NFFT__FES__ "$ & ", r); 201 202 NFFT(finalize)(&cp); 203 SOLVER(finalize_complex)(&ip); 204 NFFT(finalize)(&p); 205} 206 207/** Main routine */ 208int main(int argc, char **argv) 209{ 210 int M_cv; 211 212 if (argc < 3) 213 { 214 fprintf(stderr, "Call this program from the Matlab script glacier.m!"); 215 return EXIT_FAILURE; 216 } 217 218 if (argc == 3) 219 glacier(atoi(argv[1]), atoi(argv[2])); 220 else 221 for (M_cv = atoi(argv[3]); M_cv <= atoi(argv[5]); M_cv += atoi(argv[4])) 222 { 223 fprintf(stderr, "\nM_cv=%d,\t", M_cv); 224 printf("$%d$ & ", M_cv); 225 fprintf(stderr, "cgne+damp: "); 226 glacier_cv(atoi(argv[1]), atoi(argv[2]), M_cv, CGNE | PRECOMPUTE_DAMP); 227 //fprintf(stderr,"cgne: "); 228 //glacier_cv(atoi(argv[1]),atoi(argv[2]),M_cv,CGNE); 229 fprintf(stderr, "cgnr: "); 230 glacier_cv(atoi(argv[1]), atoi(argv[2]), M_cv, CGNR); 231 fprintf(stderr, "cgnr: "); 232 glacier_cv(atoi(argv[1]) / 4, atoi(argv[2]), M_cv, CGNR); 233 printf("XXX \\\\\n"); 234 } 235 236 fprintf(stderr, "\n"); 237 238 return EXIT_SUCCESS; 239} 240/* \} */ 241