commit last state
[ppam-mpi.git] / code / src / Algorithm / compute_coefficients.c
1 #include <stdlib.h>
2 #include "Util/types.h"
3 #include <math.h>
4 #include <gsl/gsl_wavelet.h>
5 #include <gsl/gsl_spline.h>
6
7 // compute rows of the matrix of reduced coordinates
8 void compute_coefficients(PowerCurve* powerCurves, uint32_t nbSeries, uint32_t nbValues,
9 Real* reducedCoordinates, uint32_t index, uint32_t nbReducedCoordinates)
10 {
11 uint32_t D = (1 << nbReducedCoordinates);
12 Real* x = (Real*) malloc(nbValues*sizeof(Real));
13 for (uint32_t i=0; i<nbValues; i++)
14 x[i] = i;
15 Real* y = (Real*) malloc(nbValues*sizeof(Real));
16 Real* interpolatedTranformedCurve = (Real*) malloc(D*sizeof(Real));
17 gsl_interp* linearInterpolation = gsl_interp_alloc(gsl_interp_linear, nbValues);
18 gsl_interp_accel* acc = gsl_interp_accel_alloc();
19 gsl_wavelet_workspace* work = gsl_wavelet_workspace_alloc(D);
20 //gsl_wavelet* w = gsl_wavelet_alloc(gsl_wavelet_bspline, 206); //used for starlight
21 gsl_wavelet* w = gsl_wavelet_alloc(gsl_wavelet_daubechies, 6); //used for power curves
22 //gsl_wavelet* w = gsl_wavelet_alloc(gsl_wavelet_haar, 2);
23 for (uint32_t i = 0; i < nbSeries; i++)
24 {
25 //Spline interpolation to have D = 2^u sample points
26 for (uint32_t j=0; j<nbValues; j++)
27 y[j] = powerCurves[i].values[j];
28 gsl_interp_init(linearInterpolation, x, y, nbValues);
29 for (uint32_t j=0; j<D; j++)
30 {
31 interpolatedTranformedCurve[j] =
32 gsl_interp_eval(linearInterpolation, x, y, j*((Real)(nbValues-1)/(D-1)), acc);
33 }
34 //DWT transform (in place) on interpolated curve [TODO: clarify stride parameter]
35 gsl_wavelet_transform_forward(w, interpolatedTranformedCurve, 1, D, work);
36 //Fill reducedCoordinates with energy contributions
37 uint32_t t0 = 1;
38 uint32_t t1 = 1;
39 for (uint32_t j=0; j<nbReducedCoordinates; j++)
40 {
41 t1 += (1 << j);
42 //reducedCoordinates[(index+i)*nbReducedCoordinates+j] = sqrt( sum( x[t0:t1]^2 ) )
43 Real sumOnSegment = 0.0;
44 for (uint32_t u=t0; u<t1; u++)
45 sumOnSegment += interpolatedTranformedCurve[u]*interpolatedTranformedCurve[u];
46 reducedCoordinates[(index+i)*nbReducedCoordinates+j] = sqrt(sumOnSegment);
47 t0 = t1;
48 }
49 }
50 gsl_interp_free(linearInterpolation);
51 gsl_interp_accel_free(acc);
52 free(x);
53 free(y);
54 free(interpolatedTranformedCurve);
55 gsl_wavelet_free(w);
56 gsl_wavelet_workspace_free(work);
57 }