|
@@ -31,7 +31,7 @@
|
|
|
/*-------------------------------------------------------------------------------------------*/
|
|
|
int cross_correlation(struct Map_info *Map, double passWE, double passNS)
|
|
|
/*
|
|
|
- Map: Map in which cross crorrelation will take values
|
|
|
+ Map: Vector map from which cross-crorrelation will take values
|
|
|
passWE: spline step in West-East direction
|
|
|
passNS: spline step in North-South direction
|
|
|
|
|
@@ -46,6 +46,9 @@ int cross_correlation(struct Map_info *Map, double passWE, double passNS)
|
|
|
|
|
|
/* double lambda[PARAM_LAMBDA] = { 0.0001, 0.001, 0.01, 0.1, 1.0, 10.0 }; */ /* Fixed values (by the moment) */
|
|
|
double lambda[PARAM_LAMBDA] = { 0.0001, 0.001, 0.005, 0.01, 0.02, 0.05 }; /* Fixed values (by the moment) */
|
|
|
+ /* a more exhaustive search:
|
|
|
+ #define PARAM_LAMBDA 11
|
|
|
+ double lambda[PARAM_LAMBDA] = { 0.0001, 0.0005, 0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1.0, 5.0, 10.0 }; */
|
|
|
|
|
|
double *TN, *Q, *parVect; /* Interpolation and least-square vectors */
|
|
|
double **N, **obsVect; /* Interpolation and least-square matrix */
|
|
@@ -92,7 +95,7 @@ int cross_correlation(struct Map_info *Map, double passWE, double passNS)
|
|
|
int BW;
|
|
|
double mean_reg, *obs_mean;
|
|
|
|
|
|
- int nrec, ctype = 0;
|
|
|
+ int nrec, ctype = 0, verbosity;
|
|
|
struct field_info *Fi;
|
|
|
dbDriver *driver_cats;
|
|
|
|
|
@@ -100,6 +103,8 @@ int cross_correlation(struct Map_info *Map, double passWE, double passNS)
|
|
|
rms = G_alloc_vector(PARAM_LAMBDA); /* number of parameter used used for cross validation */
|
|
|
stdev = G_alloc_vector(PARAM_LAMBDA);
|
|
|
|
|
|
+ verbosity = G_verbose(); /* store for later reset */
|
|
|
+
|
|
|
/* Working with attributes */
|
|
|
if (bspline_field > 0) {
|
|
|
db_CatValArray_init(&cvarr);
|
|
@@ -154,30 +159,31 @@ int cross_correlation(struct Map_info *Map, double passWE, double passNS)
|
|
|
TN = G_alloc_vector(nparam_spl); /* vector */
|
|
|
parVect = G_alloc_vector(nparam_spl); /* Parameters vector */
|
|
|
obsVect = G_alloc_matrix(ndata, 3); /* Observation vector */
|
|
|
- Q = G_alloc_vector(ndata); /* "a priori" var-cov matrix */
|
|
|
+ Q = G_alloc_vector(ndata); /* "a priori" var-cov matrix */
|
|
|
|
|
|
obs_mean = G_alloc_vector(ndata);
|
|
|
stat_vect = alloc_Stats(ndata);
|
|
|
|
|
|
for (lbd = 0; lbd < PARAM_LAMBDA; lbd++) { /* For each lambda value */
|
|
|
|
|
|
- G_message(_("Begining cross validation with "
|
|
|
- "lambda_i=%.4f..."), lambda[lbd]);
|
|
|
-
|
|
|
+ G_message(_("Beginning cross validation with "
|
|
|
+ "lambda_i=%.4f ... (%d of %d)"), lambda[lbd],
|
|
|
+ lbd+1, PARAM_LAMBDA);
|
|
|
+
|
|
|
/*
|
|
|
- How cross correlation algorithm is done:
|
|
|
- For each cicle, only the first ndata-1 "observ" elements are considered for the
|
|
|
+ How the cross correlation algorithm is done:
|
|
|
+ For each cycle, only the first ndata-1 "observ" elements are considered for the
|
|
|
interpolation. Within every interpolation mean is calculated to lowering edge
|
|
|
- errors. The point let out will be used for an estimation. The error between the
|
|
|
+ errors. The point left out will be used for an estimation. The error between the
|
|
|
estimation and the observation is recorded for further statistics.
|
|
|
- At the end of the cicle, the last point, that is, the ndata-1 index, and the point
|
|
|
+ At the end of the cycle, the last point, that is, the ndata-1 index, and the point
|
|
|
with j index are swapped.
|
|
|
*/
|
|
|
for (j = 0; j < ndata; j++) { /* Cross Correlation will use all ndata points */
|
|
|
- double out_x, out_y, out_z; /* This point is let out */
|
|
|
+ double out_x, out_y, out_z; /* This point is left out */
|
|
|
|
|
|
for (i = 0; i < ndata; i++) { /* Each time, only the first ndata-1 points */
|
|
|
- double dval; /* are considered in the interpolation */
|
|
|
+ double dval; /* are considered in the interpolation */
|
|
|
|
|
|
/* Setting obsVect vector & Q matrix */
|
|
|
Q[i] = 1; /* Q=I */
|
|
@@ -228,7 +234,7 @@ int cross_correlation(struct Map_info *Map, double passWE, double passNS)
|
|
|
for (i = 0; i < ndata; i++)
|
|
|
obsVect[i][2] -= mean_reg;
|
|
|
|
|
|
- /* This is let out */
|
|
|
+ /* This is left out */
|
|
|
out_x = observ[ndata - 1].coordX;
|
|
|
out_y = observ[ndata - 1].coordY;
|
|
|
out_z = obsVect[ndata - 1][2];
|
|
@@ -252,7 +258,9 @@ int cross_correlation(struct Map_info *Map, double passWE, double passNS)
|
|
|
if (bilin) interpolation (&interp, P_BILINEAR);
|
|
|
else interpolation (&interp, P_BICUBIC);
|
|
|
*/
|
|
|
+ G_set_verbose(G_verbose_min());
|
|
|
G_math_solver_cholesky_sband(N, parVect, TN, nparam_spl, BW);
|
|
|
+ G_set_verbose(verbosity);
|
|
|
|
|
|
/* Estimation of j-point */
|
|
|
if (bilin)
|
|
@@ -267,12 +275,15 @@ int cross_correlation(struct Map_info *Map, double passWE, double passNS)
|
|
|
nsplx, nsply, region.west,
|
|
|
region.south, parVect);
|
|
|
|
|
|
- /*Difference between estimated and observated i-point */
|
|
|
+ /* Difference between estimated and observated i-point */
|
|
|
stat_vect.error[j] = out_z - stat_vect.estima[j];
|
|
|
G_debug(1, "CrossCorrelation: stat_vect.error[%d] = %lf", j,
|
|
|
stat_vect.error[j]);
|
|
|
|
|
|
- observ = swap(observ, j, ndata - 1); /* Once the last value is let out, it is swap with j-value */
|
|
|
+ /* Once the last value is left out, it is swaped with j-value */
|
|
|
+ observ = swap(observ, j, ndata - 1);
|
|
|
+
|
|
|
+ G_percent(j, ndata, 2);
|
|
|
}
|
|
|
|
|
|
mean[lbd] = calc_mean(stat_vect.error, stat_vect.n_points);
|
|
@@ -282,8 +293,9 @@ int cross_correlation(struct Map_info *Map, double passWE, double passNS)
|
|
|
calc_standard_deviation(stat_vect.error, stat_vect.n_points);
|
|
|
|
|
|
G_message(_("Mean = %.5lf"), mean[lbd]);
|
|
|
- G_message(_("Root Means Square (RMS) = %.5lf"),
|
|
|
+ G_message(_("Root Mean Square (RMS) = %.5lf"),
|
|
|
rms[lbd]);
|
|
|
+ G_message("---");
|
|
|
} /* ENDFOR each lambda value */
|
|
|
|
|
|
G_free_matrix(N);
|
|
@@ -293,7 +305,7 @@ int cross_correlation(struct Map_info *Map, double passWE, double passNS)
|
|
|
G_free_vector(parVect);
|
|
|
#ifdef nodef
|
|
|
/*TODO: if the minimum lambda is wanted, the function declaration must be changed */
|
|
|
- /*At this moment, consider rms only */
|
|
|
+ /* At this moment, consider rms only */
|
|
|
rms_min = find_minimum(rms, &lbd_min);
|
|
|
stdev_min = find_minimum(stdev, &lbd_min);
|
|
|
|
|
@@ -308,11 +320,11 @@ int cross_correlation(struct Map_info *Map, double passWE, double passNS)
|
|
|
*lambda_min = lambda[lbd_min];
|
|
|
#endif
|
|
|
|
|
|
- G_message(_("Results into a table:"));
|
|
|
- G_message(_(" lambda | mean | rms |"));
|
|
|
+ G_message(_("Table of results:"));
|
|
|
+ fprintf(stdout, _(" lambda | mean | rms |\n"));
|
|
|
for (lbd = 0; lbd < PARAM_LAMBDA; lbd++) {
|
|
|
- G_message(_(" %-10.5f| %-12.4f| %-12.4f|"), lambda[lbd],
|
|
|
- mean[lbd], rms[lbd]);
|
|
|
+ fprintf(stdout, " %9.5f | %10.4f | %10.4f |\n", lambda[lbd],
|
|
|
+ mean[lbd], rms[lbd]);
|
|
|
}
|
|
|
|
|
|
G_free_vector(mean);
|
|
@@ -428,8 +440,8 @@ double find_minimum(double *values, int *l_min)
|
|
|
|
|
|
struct Point *swap(struct Point *point, int a, int b)
|
|
|
{
|
|
|
-
|
|
|
- SWAP(point[a].coordX, point[b].coordX); /* Once the last value is let out, it is swap with j-value */
|
|
|
+ /* Once the last value is left out, it is swaped with j-value */
|
|
|
+ SWAP(point[a].coordX, point[b].coordX);
|
|
|
SWAP(point[a].coordY, point[b].coordY);
|
|
|
SWAP(point[a].coordZ, point[b].coordZ);
|
|
|
SWAP(point[a].cat, point[b].cat);
|