iterS2.c
/*
* Copyright (c) 2002, 2017 Jens Keiner, Stefan Kunis, Daniel Potts
*
* This program is free software; you can redistribute it and/or modify it under
* the terms of the GNU General Public License as published by the Free Software
* Foundation; either version 2 of the License, or (at your option) any later
* version.
*
* This program is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
* FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
* details.
*
* You should have received a copy of the GNU General Public License along with
* this program; if not, write to the Free Software Foundation, Inc., 51
* Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
/* iterS2 - Iterative reconstruction on the sphere S2 */
#include "config.h"
/* Include standard C headers. */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#ifdef HAVE_COMPLEX_H
#include <complex.h>
#endif
/* Include NFFT 3 utilities headers. */
/* Include NFFT3 library header. */
#include "nfft3.h"
#include "infft.h"
#include "legendre.h"
static void voronoi_weights_S2(R *w, R *xi, INT M)
{
R *x;
R *y;
R *z;
int j;
int k;
int el;
int Mlocal = M;
int lnew;
int ier;
int *list;
int *lptr;
int *lend;
int *near;
int *next;
R *dist;
int *ltri;
int *listc;
int nb;
R *xc;
R *yc;
R *zc;
R *rc;
R *vr;
int lp;
int lpl;
int kv;
R a;
/* Allocate memory for auxilliary arrays. */
x = (R*)X(malloc)(M * sizeof(R));
y = (R*)X(malloc)(M * sizeof(R));
z = (R*)X(malloc)(M * sizeof(R));
list = (int*)X(malloc)((6*M-12+1)*sizeof(int));
lptr = (int*)X(malloc)((6*M-12+1)*sizeof(int));
lend = (int*)X(malloc)((M+1)*sizeof(int));
near = (int*)X(malloc)((M+1)*sizeof(int));
next = (int*)X(malloc)((M+1)*sizeof(int));
dist = (R*)X(malloc)((M+1)*sizeof(R));
ltri = (int*)X(malloc)((6*M+1)*sizeof(int));
listc = (int*)X(malloc)((6*M-12+1)*sizeof(int));
xc = (R*)X(malloc)((2*M-4+1)*sizeof(R));
yc = (R*)X(malloc)((2*M-4+1)*sizeof(R));
zc = (R*)X(malloc)((2*M-4+1)*sizeof(R));
rc = (R*)X(malloc)((2*M-4+1)*sizeof(R));
vr = (R*)X(malloc)(3*(2*M-4+1)*sizeof(R));
/* Convert from spherical Coordinates in [0,1/2]x[-1/2,1/2) to Cartesian
* coordinates. */
for (k = 0; k < M; k++)
{
x[k] = SIN(K2PI*xi[2*k+1])*COS(K2PI*xi[2*k]);
y[k] = SIN(K2PI*xi[2*k+1])*SIN(K2PI*xi[2*k]);
z[k] = COS(K2PI*xi[2*k+1]);
}
/* Generate Delaunay triangulation. */
trmesh_(&Mlocal, x, y, z, list, lptr, lend, &lnew, near, next, dist, &ier);
/* Check error flag. */
if (ier == 0)
{
/* Get Voronoi vertices. */
crlist_(&Mlocal, &Mlocal, x, y, z, list, lend, lptr, &lnew, ltri, listc, &nb, xc,
yc, zc, rc, &ier);
if (ier == 0)
{
/* Calcuate sizes of Voronoi regions. */
for (k = 0; k < M; k++)
{
/* Get last neighbour index. */
lpl = lend[k];
lp = lpl;
j = 0;
vr[3*j] = x[k];
vr[3*j+1] = y[k];
vr[3*j+2] = z[k];
do
{
j++;
/* Get next neighbour. */
lp = lptr[lp-1];
kv = listc[lp-1];
vr[3*j] = xc[kv-1];
vr[3*j+1] = yc[kv-1];
vr[3*j+2] = zc[kv-1];
/* fprintf(stderr, "lp = %ld\t", lp); */
} while (lp != lpl);
a = 0;
for (el = 0; el < j; el++)
{
a += areas_(vr, &vr[3*(el+1)],&vr[3*(((el+1)%j)+1)]);
}
w[k] = a;
}
}
}
/* Deallocate memory. */
X(free)(x);
X(free)(y);
X(free)(z);
X(free)(list);
X(free)(lptr);
X(free)(lend);
X(free)(near);
X(free)(next);
X(free)(dist);
X(free)(ltri);
X(free)(listc);
X(free)(xc);
X(free)(yc);
X(free)(zc);
X(free)(rc);
X(free)(vr);
}
enum boolean {NO = 0, YES = 1};
int main (int argc, char **argv)
{
int T;
int N;
int M;
int M2;
int t; /* Index variable for testcases */
nfsft_plan plan; /* NFSFT plan */
nfsft_plan plan2; /* NFSFT plan */
solver_plan_complex iplan; /* NFSFT plan */
int j; /* */
int k; /* */
int m; /* */
int use_nfsft; /* */
int use_nfft; /* */
int use_fpt; /* */
int cutoff;
double threshold;
double re;
double im;
double a;
double xs;
double *ys;
double *temp;
double _Complex *temp2;
int qlength;
double *qweights;
fftw_plan fplan;
fpt_set set;
int npt;
int npt_exp;
double *alpha, *beta, *gamma;
/* Read the number of testcases. */
fscanf(stdin,"testcases=%d\n",&T);
fprintf(stderr,"%d\n",T);
/* Process each testcase. */
for (t = 0; t < T; t++)
{
/* Check if the fast transform shall be used. */
fscanf(stdin,"nfsft=%d\n",&use_nfsft);
fprintf(stderr,"%d\n",use_nfsft);
if (use_nfsft != NO)
{
/* Check if the NFFT shall be used. */
fscanf(stdin,"nfft=%d\n",&use_nfft);
fprintf(stderr,"%d\n",use_nfsft);
if (use_nfft != NO)
{
/* Read the cut-off parameter. */
fscanf(stdin,"cutoff=%d\n",&cutoff);
fprintf(stderr,"%d\n",cutoff);
}
else
{
/* TODO remove this */
/* Initialize unused variable with dummy value. */
cutoff = 1;
}
/* Check if the fast polynomial transform shall be used. */
fscanf(stdin,"fpt=%d\n",&use_fpt);
fprintf(stderr,"%d\n",use_fpt);
if (use_fpt != NO)
{
/* Read the NFSFT threshold parameter. */
fscanf(stdin,"threshold=%lf\n",&threshold);
fprintf(stderr,"%lf\n",threshold);
}
else
{
/* TODO remove this */
/* Initialize unused variable with dummy value. */
threshold = 1000.0;
}
}
else
{
/* TODO remove this */
/* Set dummy values. */
use_nfft = NO;
use_fpt = NO;
cutoff = 3;
threshold = 1000.0;
}
/* Read the bandwidth. */
fscanf(stdin,"bandwidth=%d\n",&N);
fprintf(stderr,"%d\n",N);
/* Do precomputation. */
nfsft_precompute(N,threshold,
((use_nfsft==NO)?(NFSFT_NO_FAST_ALGORITHM):(0U/*NFSFT_NO_DIRECT_ALGORITHM*/)), 0U);
/* Read the number of nodes. */
fscanf(stdin,"nodes=%d\n",&M);
fprintf(stderr,"%d\n",M);
/* */
if ((N+1)*(N+1) > M)
{
X(next_power_of_2_exp_int)(N, &npt, &npt_exp);
fprintf(stderr, "npt = %d, npt_exp = %d\n", npt, npt_exp);
fprintf(stderr,"Optimal interpolation!\n");
ys = (double*) nfft_malloc((N+1)*sizeof(double));
temp = (double*) nfft_malloc((2*N+1)*sizeof(double));
temp2 = (double _Complex*) nfft_malloc((N+1)*sizeof(double _Complex));
a = 0.0;
for (j = 0; j <= N; j++)
{
xs = 2.0 + (2.0*j)/(N+1);
ys[j] = (2.0-((j == 0)?(1.0):(0.0)))*4.0*nfft_bsplines(4,xs);
//fprintf(stdout,"%3d: g(%le) = %le\n",j,xs,ys[j]);
a += ys[j];
}
//fprintf(stdout,"a = %le\n",a);
for (j = 0; j <= N; j++)
{
ys[j] *= 1.0/a;
}
qlength = 2*N+1;
qweights = (double*) nfft_malloc(qlength*sizeof(double));
fplan = fftw_plan_r2r_1d(N+1, qweights, qweights, FFTW_REDFT00, 0U);
for (j = 0; j < N+1; j++)
{
qweights[j] = -2.0/(4*j*j-1);
}
fftw_execute(fplan);
qweights[0] *= 0.5;
for (j = 0; j < N+1; j++)
{
qweights[j] *= 1.0/(2.0*N+1.0);
qweights[2*N+1-1-j] = qweights[j];
}
fplan = fftw_plan_r2r_1d(2*N+1, temp, temp, FFTW_REDFT00, 0U);
for (j = 0; j <= N; j++)
{
temp[j] = ((j==0 || j == 2*N)?(1.0):(0.5))*ys[j];
}
for (j = N+1; j < 2*N+1; j++)
{
temp[j] = 0.0;
}
fftw_execute(fplan);
for (j = 0; j < 2*N+1; j++)
{
temp[j] *= qweights[j];
}
fftw_execute(fplan);
for (j = 0; j < 2*N+1; j++)
{
temp[j] *= ((j==0 || j == 2*N)?(1.0):(0.5));
if (j <= N)
{
temp2[j] = temp[j];
}
}
set = fpt_init(1, npt_exp, 0U);
alpha = (double*) nfft_malloc((N+2)*sizeof(double));
beta = (double*) nfft_malloc((N+2)*sizeof(double));
gamma = (double*) nfft_malloc((N+2)*sizeof(double));
alpha_al_row(alpha, N, 0);
beta_al_row(beta, N, 0);
gamma_al_row(gamma, N, 0);
fpt_precompute(set, 0, alpha, beta, gamma, 0, 1000.0);
fpt_transposed(set,0, temp2, temp2, N, 0U);
fpt_finalize(set);
nfft_free(alpha);
nfft_free(beta);
nfft_free(gamma);
fftw_destroy_plan(fplan);
nfft_free(qweights);
nfft_free(ys);
nfft_free(temp);
}
/* Init transform plans. */
nfsft_init_guru(&plan, N, M,
((use_nfft!=0)?(0U):(NFSFT_USE_NDFT)) |
((use_fpt!=0)?(0U):(NFSFT_USE_DPT)) | NFSFT_MALLOC_F | NFSFT_MALLOC_X |
NFSFT_MALLOC_F_HAT | NFSFT_NORMALIZED | NFSFT_ZERO_F_HAT,
PRE_PHI_HUT | PRE_PSI | FFTW_INIT |
FFT_OUT_OF_PLACE,
cutoff);
if ((N+1)*(N+1) > M)
{
solver_init_advanced_complex(&iplan, (nfft_mv_plan_complex*)(&plan), CGNE | PRECOMPUTE_DAMP);
}
else
{
solver_init_advanced_complex(&iplan, (nfft_mv_plan_complex*)(&plan), CGNR | PRECOMPUTE_WEIGHT | PRECOMPUTE_DAMP);
}
/* Read the nodes and function values. */
for (j = 0; j < M; j++)
{
fscanf(stdin,"%le %le %le %le\n",&plan.x[2*j+1],&plan.x[2*j],&re,&im);
plan.x[2*j+1] = plan.x[2*j+1]/(2.0*KPI);
plan.x[2*j] = plan.x[2*j]/(2.0*KPI);
if (plan.x[2*j] >= 0.5)
{
plan.x[2*j] = plan.x[2*j] - 1;
}
iplan.y[j] = re + _Complex_I * im;
fprintf(stderr,"%le %le %le %le\n",plan.x[2*j+1],plan.x[2*j],
creal(iplan.y[j]),cimag(iplan.y[j]));
}
/* Read the number of nodes. */
fscanf(stdin,"nodes_eval=%d\n",&M2);
fprintf(stderr,"%d\n",M2);
/* Init transform plans. */
nfsft_init_guru(&plan2, N, M2,
((use_nfft!=0)?(0U):(NFSFT_USE_NDFT)) |
((use_fpt!=0)?(0U):(NFSFT_USE_DPT)) | NFSFT_MALLOC_F | NFSFT_MALLOC_X |
NFSFT_MALLOC_F_HAT | NFSFT_NORMALIZED | NFSFT_ZERO_F_HAT,
PRE_PHI_HUT | PRE_PSI | FFTW_INIT |
FFT_OUT_OF_PLACE,
cutoff);
/* Read the nodes and function values. */
for (j = 0; j < M2; j++)
{
fscanf(stdin,"%le %le\n",&plan2.x[2*j+1],&plan2.x[2*j]);
plan2.x[2*j+1] = plan2.x[2*j+1]/(2.0*KPI);
plan2.x[2*j] = plan2.x[2*j]/(2.0*KPI);
if (plan2.x[2*j] >= 0.5)
{
plan2.x[2*j] = plan2.x[2*j] - 1;
}
fprintf(stderr,"%le %le\n",plan2.x[2*j+1],plan2.x[2*j]);
}
nfsft_precompute_x(&plan);
nfsft_precompute_x(&plan2);
/* Frequency weights. */
if ((N+1)*(N+1) > M)
{
/* Compute Voronoi weights. */
//voronoi_weights_S2(iplan.w, plan.x, M);
/* Print out Voronoi weights. */
/*a = 0.0;
for (j = 0; j < plan.M_total; j++)
{
fprintf(stderr,"%le\n",iplan.w[j]);
a += iplan.w[j];
}
fprintf(stderr,"sum = %le\n",a);*/
for (j = 0; j < plan.N_total; j++)
{
iplan.w_hat[j] = 0.0;
}
for (k = 0; k <= N; k++)
{
for (j = -k; j <= k; j++)
{
iplan.w_hat[NFSFT_INDEX(k,j,&plan)] = 1.0/(pow(k+1.0,2.0)); /*temp2[j]*/;
}
}
}
else
{
for (j = 0; j < plan.N_total; j++)
{
iplan.w_hat[j] = 0.0;
}
for (k = 0; k <= N; k++)
{
for (j = -k; j <= k; j++)
{
iplan.w_hat[NFSFT_INDEX(k,j,&plan)] = 1/(pow(k+1.0,2.5));
}
}
/* Compute Voronoi weights. */
voronoi_weights_S2(iplan.w, plan.x, M);
/* Print out Voronoi weights. */
a = 0.0;
for (j = 0; j < plan.M_total; j++)
{
fprintf(stderr,"%le\n",iplan.w[j]);
a += iplan.w[j];
}
fprintf(stderr,"sum = %le\n",a);
}
fprintf(stderr, "N_total = %d\n", plan.N_total);
fprintf(stderr, "M_total = %d\n", plan.M_total);
/* init some guess */
for (k = 0; k < plan.N_total; k++)
{
iplan.f_hat_iter[k] = 0.0;
}
/* inverse trafo */
solver_before_loop_complex(&iplan);
/*for (k = 0; k < plan.M_total; k++)
{
printf("%le %le\n",creal(iplan.r_iter[k]),cimag(iplan.r_iter[k]));
}*/
for (m = 0; m < 29; m++)
{
fprintf(stderr,"Residual ||r||=%e,\n",sqrt(iplan.dot_r_iter));
solver_loop_one_step_complex(&iplan);
}
/*CSWAP(iplan.f_hat_iter, plan.f_hat);
nfsft_trafo(&plan);
CSWAP(iplan.f_hat_iter, plan.f_hat);
a = 0.0;
b = 0.0;
for (k = 0; k < plan.M_total; k++)
{
printf("%le %le %le\n",cabs(iplan.y[k]),cabs(plan.f[k]),
cabs(iplan.y[k]-plan.f[k]));
a += cabs(iplan.y[k]-plan.f[k])*cabs(iplan.y[k]-plan.f[k]);
b += cabs(iplan.y[k])*cabs(iplan.y[k]);
}
fprintf(stderr,"relative error in 2-norm: %le\n",a/b);*/
CSWAP(iplan.f_hat_iter, plan2.f_hat);
nfsft_trafo(&plan2);
CSWAP(iplan.f_hat_iter, plan2.f_hat);
for (k = 0; k < plan2.M_total; k++)
{
fprintf(stdout,"%le\n",cabs(plan2.f[k]));
}
solver_finalize_complex(&iplan);
nfsft_finalize(&plan);
nfsft_finalize(&plan2);
/* Delete precomputed data. */
nfsft_forget();
if ((N+1)*(N+1) > M)
{
nfft_free(temp2);
}
} /* Process each testcase. */
/* Return exit code for successful run. */
return EXIT_SUCCESS;
}
/* \} */
X
#define X(name)
Include header for C99 complex datatype.
Definition fastsum.h:57
main
int main(int argc, char **argv)
The main program.
Definition iterS2.c:181
fpt_init
fpt_set fpt_init(const int M, const int t, const unsigned int flags)
Definition fpt.c:795
fpt_precompute
void fpt_precompute(fpt_set set, const int m, double *alpha, double *beta, double *gam, int k_start, const double threshold)
Definition fpt.c:1307
FFT_OUT_OF_PLACE
#define FFT_OUT_OF_PLACE
Definition nfft3.h:190
PRE_PSI
#define PRE_PSI
Definition nfft3.h:185
FFTW_INIT
#define FFTW_INIT
Definition nfft3.h:191
PRE_PHI_HUT
#define PRE_PHI_HUT
Definition nfft3.h:181
CSWAP
#define CSWAP(x, y)
Swap two vectors.
Definition infft.h:139
NFSFT_MALLOC_X
#define NFSFT_MALLOC_X
Definition nfft3.h:582
nfsft_forget
void nfsft_forget(void)
Definition nfsft.c:579
nfsft_trafo
void nfsft_trafo(nfsft_plan *plan)
Definition nfsft.c:1068
NFSFT_NORMALIZED
#define NFSFT_NORMALIZED
Definition nfft3.h:579
NFSFT_USE_DPT
#define NFSFT_USE_DPT
Definition nfft3.h:581
NFSFT_ZERO_F_HAT
#define NFSFT_ZERO_F_HAT
Definition nfft3.h:596
NFSFT_INDEX
#define NFSFT_INDEX(k, n, plan)
Definition nfft3.h:599
NFSFT_NO_FAST_ALGORITHM
#define NFSFT_NO_FAST_ALGORITHM
Definition nfft3.h:595
nfsft_finalize
void nfsft_finalize(nfsft_plan *plan)
Definition nfsft.c:625
NFSFT_MALLOC_F_HAT
#define NFSFT_MALLOC_F_HAT
Definition nfft3.h:583
NFSFT_USE_NDFT
#define NFSFT_USE_NDFT
Definition nfft3.h:580
nfsft_precompute
void nfsft_precompute(int N, double kappa, unsigned int nfsft_flags, unsigned int fpt_flags)
Definition nfsft.c:376
NFSFT_MALLOC_F
#define NFSFT_MALLOC_F
Definition nfft3.h:584
CGNR
#define CGNR
Definition nfft3.h:802
PRECOMPUTE_DAMP
#define PRECOMPUTE_DAMP
Definition nfft3.h:806
PRECOMPUTE_WEIGHT
#define PRECOMPUTE_WEIGHT
Definition nfft3.h:805
CGNE
#define CGNE
Definition nfft3.h:803
infft.h
Internal header file for auxiliary definitions and functions.
nfft3.h
Header file for the nfft3 library.
fpt_set
struct fpt_set_s_ * fpt_set
A set of precomputed data for a set of DPT transforms of equal maximum length.
Definition nfft3.h:633
nfft_malloc
void * nfft_malloc(size_t n)
nfft_free
void nfft_free(void *p)
nfft_mv_plan_complex
Definition nfft3.h:177
nfsft_plan
data structure for an NFSFT (nonequispaced fast spherical Fourier transform) plan with double precisi...
Definition nfft3.h:575
nfsft_plan::N_total
NFFT_INT N_total
Total number of Fourier coefficients.
Definition nfft3.h:575
nfsft_plan::M_total
NFFT_INT M_total
Total number of samples.
Definition nfft3.h:575
nfsft_plan::f
fftw_complex * f
Samples.
Definition nfft3.h:575
nfsft_plan::x
double * x
the nodes for ,
Definition nfft3.h:575
nfsft_plan::f_hat
fftw_complex * f_hat
Fourier coefficients.
Definition nfft3.h:575
solver_plan_complex
data structure for an inverse NFFT plan with double precision
Definition nfft3.h:796
solver_plan_complex::w
double * w
weighting factors
Definition nfft3.h:796
solver_plan_complex::w_hat
double * w_hat
damping factors
Definition nfft3.h:796
solver_plan_complex::dot_r_iter
double dot_r_iter
weighted dotproduct of r_iter
Definition nfft3.h:796
solver_plan_complex::y
fftw_complex * y
right hand side, samples
Definition nfft3.h:796
solver_plan_complex::f_hat_iter
fftw_complex * f_hat_iter
iterative solution
Definition nfft3.h:796