https://github.com/cran/gstat
Tip revision: 22b93811cf4a7d81cd4071888f77de9088d8bdb9 authored by Edzer J. Pebesma on 12 November 2008, 16:36:19 UTC
version 0.9-54
version 0.9-54
Tip revision: 22b9381
getest.c
/*
Gstat, a program for geostatistical modelling, prediction and simulation
Copyright 1992, 2003 (C) Edzer J. Pebesma
Edzer J. Pebesma, e.pebesma@geo.uu.nl
Department of physical geography, Utrecht University
P.O. Box 80.115, 3508 TC Utrecht, The Netherlands
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. As a special exception, linking
this program with the Qt library is permitted.
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., 675 Mass Ave, Cambridge, MA 02139, USA.
(read also the files COPYING and Copyright)
*/
/*
* getest.c: choose the predicion function at one prediction location
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <limits.h>
#include "defs.h"
#include "data.h"
#include "utils.h"
#include "mapio.h"
#include "userio.h"
#include "debug.h"
#include "vario.h"
#include "sem.h"
#include "fit.h"
#include "glvars.h"
#include "lm.h"
#include "gls.h"
#include "sim.h"
#include "msim.h"
#include "stat.h"
#include "block.h"
#include "polygon.h"
#include "getest.h"
static void est_quantile_div(DATA *data, double *est, int div);
static void est_skew_kurt(DATA *data, double *est);
static double inverse_dist(DATA *data, DPOINT *where, double idPow);
static void save_variogram_parameters(VARIOGRAM *v);
static void reset_variogram_parameters(VARIOGRAM *v);
static double *vgm_pars = NULL;
static void calc_poly(DPOINT *where, double *est);
void get_est(DATA **data, METHOD method, DPOINT *where, double *est) {
/*
* given all data[i]->sel, get_est returns one or more values to *est,
* according to the method of calculation set in *method
* sel must be not NULL, where contains the location of the estimate
* PRE: data, where, est;
* USES: several estimation routines.
*/
DPOINT *block = NULL;
VARIOGRAM *v;
int i, j, n_vars, n_sel, *is_pt;
static GRIDMAP **m = NULL;
unsigned int row = 0, col = 0;
double *X_ori = NULL, *local_sim;
enum GLS_WHAT gls_mode = GLS_BLUP;
const double *sim = NULL; /* return value of cond_sim() */
for (i = 0; i < get_n_outfile(); i++)
set_mv_double(&est[i]);
block = get_block_p();
if (get_mode() == NSP)
ErrMsg(ER_IMPOSVAL, "Getest(): mode not specified");
if (block->x > 0.0 || block->y > 0.0 || block->z > 0.0 || get_data_area())
SET_BLOCK(where);
else
SET_POINT(where);
n_vars = get_n_vars();
if (get_mode() == STRATIFY &&
(where->u.stratum < 0 || where->u.stratum >= n_vars))
return;
local_sim = (double *) emalloc(n_vars * sizeof(double));
is_pt = (int *) emalloc(n_vars * sizeof(int));
for (i = 0; i < n_vars; i++) {
set_mv_double(&(local_sim[i]));
is_pt[i] = 0;
}
if (DEBUG_COV) {
printlog("we're at location X: %g Y: %g Z: %g\n",
where->x, where->y, where->z);
if (IS_BLOCK(where)) {
if (get_data_area())
printlog("block set in area()\n");
else
printlog("block size: dx: %g dy: %g dz: %g\n",
block->x, block->y, block->z);
} else
printlog("zero block size\n");
if (get_mode() == STRATIFY)
printlog("stratum: %d\n", where->u.stratum);
}
switch (method) {
case DIV: /* BREAKTHROUGH */
case MED:
if (get_mode() == STRATIFY)
est_quantile_div(data[where->u.stratum], est, method == DIV);
else
for (i = 0; i < n_vars; i++)
est_quantile_div(data[i], &(est[2*i]), method == DIV);
break;
case SKEW:
if (get_mode() == STRATIFY)
est_skew_kurt(data[where->u.stratum], est);
else
for (i = 0; i < n_vars; i++)
est_skew_kurt(data[i], &(est[2*i]));
break;
case IDW:
if (gl_idp < 0.0)
ErrMsg(ER_RANGE, "idp must be non-negative");
if (get_mode() == STRATIFY) {
if (data[where->u.stratum]->n_sel > 0)
est[0] = inverse_dist(data[where->u.stratum], where, gl_idp);
} else {
for (i = 0; i < n_vars; i++)
if (data[i]->n_sel > 0)
est[2 * i] = inverse_dist(data[i], where, gl_idp);
}
break;
case LSLM:
if (n_variograms_set()) {
switch(get_mode()) {
case SIMPLE:
X_ori = where->X; /* remember... */
for (i = 0; i < n_vars; i++) {
if (data[i]->n_sel > 0)
gls(&data[i], 1, GLS_BLUE, where,
&(est[2 * i]));
where->X += data[i]->n_X;
}
where->X = X_ori; /* and put back */
break;
case STRATIFY:
X_ori = where->X; /* remember... */
for (i = 0; i < where->u.stratum; i++)
where->X += data[i]->n_X;
if (data[where->u.stratum]->n_sel > 0) {
where->X += where->u.stratum; /* strictly unnecessary */
gls(&data[where->u.stratum], 1, GLS_BLUE, where, est);
}
where->X = X_ori; /* and put back */
break;
case MULTIVARIABLE:
gls(data, n_vars, GLS_BLUE, where, est);
break;
default: ErrMsg(ER_IMPOSVAL, "Getest(): wrong mode");
break;
}
} else {
switch(get_mode()) {
case STRATIFY:
X_ori = where->X; /* remember... */
for (i = 0; i < where->u.stratum; i++)
where->X += data[i]->n_X;
if (data[where->u.stratum]->n_sel > 0)
pred_lm(&data[where->u.stratum], 1, where, est);
where->X = X_ori; /* and put back */
break;
case MULTIVARIABLE:
pred_lm(data, n_vars, where, est);
break;
case SIMPLE:
X_ori = where->X; /* remember... */
for (i = 0; i < n_vars; i++) {
if (data[i]->n_sel > 0) {
pred_lm(&data[i], 1, where, &(est[2*i]));
where->X += data[i]->n_X;
}
}
where->X = X_ori; /* and put back */
break;
default: ErrMsg(ER_IMPOSVAL, "Getest(): wrong mode");
break;
}
} /* else */
break;
case SKR: /* BREAK THROUGH (SK handled in gls()) */
case OKR: /* BREAK THROUGH: (OK is special case of UK) */
case UKR:
if (method == SKR)
gls_mode = GLS_BLP;
else
gls_mode = GLS_BLUP;
switch (get_mode()) {
case SIMPLE:
X_ori = where->X; /* remember... */
for (i = 0; i < n_vars; i++) {
gls(&data[i], 1, gls_mode, where, &(est[2 * i]));
where->X += data[i]->n_X;
}
where->X = X_ori; /* and put back */
break;
case STRATIFY:
X_ori = where->X; /* remember... */
for (i = 0; i < where->u.stratum; i++)
where->X += data[i]->n_X;
gls(&data[where->u.stratum], 1, gls_mode, where, est);
where->X = X_ori; /* and put back */
break;
case MULTIVARIABLE:
gls(data, n_vars, gls_mode, where, est);
break;
default: ErrMsg(ER_IMPOSVAL, "Getest(): wrong mode"); break;
}
if (gl_order) /* order relation violations: */
correct_orv(est, n_vars, gl_order);
break;
case GSI: /* BREAK TRHOUGH: */
case ISI:
switch (get_mode()) {
case SIMPLE: /* estimates go to est[0+1],est[2+3],... */
X_ori = where->X; /* remember... */
for (i = 0; i < gl_nsim; i++) {
restore_data_sel(data, i, n_vars);
set_beta(data, i, n_vars, method);
for (j = 0; j < n_vars; j++) {
/* estimate or update: */
if (i == 0) /* estimate */
gls(&data[j], 1, data[j]->n_sel < gl_n_uk ?
GLS_BLP : GLS_BLUP, where, &(est[2*j]));
else /* update */
gls(&data[j], 1, UPDATE, where, &(est[2*j]));
/* new simulation: */
if (gl_order <= 1) {
sim = cond_sim(&(est[2*j]), 1, method,
&(is_pt[j]), gl_order);
/* put data from simulation i; */
local_sim[j] = sim[0];
}
where->X += data[j]->n_X;
}
where->X = X_ori; /* reset */
if (gl_order > 1) { /* 2,3 if pdf, 4 if cdf: */
sim = cond_sim(est, n_vars, method, is_pt, gl_order);
for (j = 0; j < n_vars; j++)
local_sim[j] = sim[j];
}
save_sim(data, where, i, n_vars, local_sim, is_pt);
}
X_ori = where->X; /* remember... */
for (i = 0; i < n_vars; i++) {
if (! is_pt[i]) {
where->attr = est[2 * i] = local_sim[i];
/* last simulation */
push_point(data[i], where);
} else
data[i]->nsim_at_data++;
where->X += data[i]->n_X;
}
where->X = X_ori; /* reset */
break;
case STRATIFY: /* estimate goes to est[0] and est[1] */
X_ori = where->X; /* remember... */
for (i = 0; i < where->u.stratum; i++)
where->X += data[i]->n_X;
for (i = 0; i < gl_nsim; i++) {
/* load data from simulation i; */
restore_data_sel(&data[where->u.stratum], i, 0);
set_beta(&data[where->u.stratum], i, 0, method);
if (i == 0) /* estimate: */
gls(&data[where->u.stratum], 1, data[where->u.stratum]->n_sel <
gl_n_uk ? GLS_BLP : GLS_BLUP, where, est);
else /* update */
gls(&data[where->u.stratum], 1, UPDATE, where, est);
sim = cond_sim(est, 1, method, is_pt, gl_order);
/* save data from simulation i; */
save_sim_strat(data[where->u.stratum], where, i,
sim[0], is_pt[0]);
}
/* store & push */
if (! is_pt[0]) {
where->attr = est[0] = sim[0]; /* last simulation */
push_point(data[where->u.stratum], where);
} else
data[where->u.stratum]->nsim_at_data++;
where->X = X_ori; /* and put back */
break;
case MULTIVARIABLE:
for (i = 0, n_sel = 0; i < n_vars; i++)
n_sel += data[i]->n_sel;
if (n_sel < gl_n_uk)
gls_mode = GLS_BLP;
else
gls_mode = GLS_BLUP;
for (i = 0; i < gl_nsim; i++) {
/* load data from simulation i; */
restore_data_sel(data, i, n_vars);
set_beta(data, i, n_vars, method);
/* estimate or update: */
if (i == 0)
gls(data, n_vars, gls_mode, where, est);
else
gls(data, n_vars, UPDATE, where, est);
/* new simulation: */
sim = cond_sim(est, n_vars, method, is_pt, gl_order);
/* put data from simulation i; */
save_sim(data, where, i, n_vars, sim, is_pt);
}
X_ori = where->X; /* remember... */
for (i = 0; i < n_vars; i++) {
if (! is_pt[i]) {
est[2 * i] = where->attr = sim[i]; /* last simulation */
push_point(data[i], where);
} else
data[i]->nsim_at_data++;
where->X += data[i]->n_X;
}
where->X = X_ori; /* and put back */
break;
default: ErrMsg(ER_IMPOSVAL, "Getest(): wrong mode"); break;
} /* switch (get_mode()) */
for (i = 0; i < n_vars; i++)
set_mv_double(&(est[2 * i + 1]));
for (i = 2 * n_vars; i < get_n_outfile(); i++)
set_mv_double(&(est[i]));
/* print_sim(); */
break;
case MAPVALUE:
if (get_n_masks() == 0)
ErrMsg(ER_VARNOTSET, "define at least one mask");
if (m == NULL) {
if (get_n_masks() > get_n_outfile())
pr_warning("define more (dummy) points(..) to get all masks");
m = (GRIDMAP **) emalloc(get_n_masks() * sizeof(GRIDMAP *));
for (i = 0; i < get_n_masks(); i++) {
m[i] = new_map(READ_ONLY);
m[i]->filename = string_dup(get_mask_name(i));
if (map_read(m[i]) == NULL)
ErrMsg(ER_READ, get_mask_name(i));
}
}
if (map_xy2rowcol(m[0], where->x, where->y, &row, &col) == 0)
for (i = 0; i < get_n_masks() && i < get_n_outfile(); i++) {
if (!map_cell_is_mv(m[i], row, col))
est[i] = map_get_cell(m[i], row, col);
}
break;
case NRS:
if (get_mode() != STRATIFY) {
for (i = 0; i < n_vars; i++) {
est[2 * i] = (double) data[i]->n_sel;
est[2 * i + 1] = (double) data[i]->oct_filled;
}
} else {
est[0] = (double) data[where->u.stratum]->n_sel;
est[1] = (double) data[where->u.stratum]->oct_filled;
}
break;
case XYP:
est[0] = where->x;
est[1] = where->y;
break;
case POLY:
calc_poly(where, est);
break;
case SPREAD:
if (get_mode() == STRATIFY) {
if (data[where->u.stratum]->n_sel > 0) {
if (data[where->u.stratum]->what_is_u != U_ISDIST)
ErrMsg(ER_IMPOSVAL, "getest() needs distances");
n_sel = data[where->u.stratum]->n_sel;
est[0] = sqrt(data[where->u.stratum]->sel[0]->u.dist2);
est[1] = sqrt(data[where->u.stratum]->sel[n_sel - 1]->u.dist2);
}
} else {
for (i = 0; i < n_vars; i++) {
if (data[i]->n_sel > 0) {
if (data[i]->what_is_u != U_ISDIST)
ErrMsg(ER_IMPOSVAL, "getest() needs distances");
n_sel = data[i]->n_sel;
est[2 * i] = sqrt(data[i]->sel[0]->u.dist2);
est[2 * i + 1] = sqrt(data[i]->sel[n_sel - 1]->u.dist2);
}
}
}
break;
case LSEM:
push_gstat_progress_handler(no_progress);
v = get_vgm(LTI(0,0));
assert(v);
v->id1 = v->id2 = v->id = 0;
v->ev->evt = SEMIVARIOGRAM;
if (data[0]->is_residual == 0)
data[0]->is_residual = 1; /* !! */
v->ev->recalc = 1;
calc_variogram(v, NULL);
if (gl_fit) {
v->ev->refit = 1;
v->ev->fit = gl_fit;
save_variogram_parameters(v);
fit_variogram(v);
/* write back locally fitted variogram model parameters: */
for (i = j = 0; i < MIN(get_n_vars(), v->n_models); i++) {
est[2 * i + 0] = v->part[i].sill;
est[2 * i + 1] = v->part[i].range[0];
j += 2;
}
if (j < get_n_vars())
est[j] = 1.0 * v->fit_is_singular;
reset_variogram_parameters(v);
update_variogram(v); /* not sure if this is needed */
} else {
/* write back locally estimated sample variogram values: */
for (i = 0; i < MIN(get_n_vars(), v->ev->n_est); i++) {
est[2 * i + 0] = v->ev->gamma[i];
est[2 * i + 1] = 1.0 * v->ev->nh[i];
}
}
pop_gstat_progress_handler();
break;
case NSP: /* FALLTRHOUGH: */
default:
ErrMsg(ER_IMPOSVAL, "getest(): method not specified");
} /* switch */
efree(local_sim);
efree(is_pt);
return;
}
static void est_quantile_div(DATA *data, double *est, int div) {
static double *list = NULL;
static int i, size = 0;
double mod = -9999.0;
int n, run, longest_run = 0;
if (data->n_sel > size)
list = (double *) erealloc(list, (size = data->n_sel) * sizeof(double));
for (i = 0; i < data->n_sel; i++)
list[i] = data->sel[i]->attr;
qsort(list, (size_t) data->n_sel, sizeof(double),
(int CDECL (*)(const void *, const void *)) d_cmp);
if (div) { /* get diversity and modus in sel: */
n = data->n_sel;
for (i = 0; i < data->n_sel - 2; i += run) {
run = 1;
while (i + run < data->n_sel && list[i] == list[i + run]) {
run++;
n--;
}
if (run > longest_run) {
longest_run = run;
mod = list[i];
}
}
est[0] = 1.0 * n;
est[1] = mod;
} else {
if (data->n_sel < 2)
return;
if (gl_quantile == 0.5)
est[0] = est[1] = est_quant(list, gl_quantile, data->n_sel);
else {
est[0] = est_quant(list, gl_quantile, data->n_sel);
est[1] = est_quant(list, 1.0 - gl_quantile, data->n_sel);
}
}
return;
}
static void est_skew_kurt(DATA *data, double *est) {
static double *list = NULL;
double mean, std, skewness = 0.0, kurtosis = 0.0, d;
static int i, size = 0;
if (data->n_sel <= 1) /* need at least 2 for variance */
return;
if (data->n_sel > size)
list = (double *) erealloc(list, (size = data->n_sel) * sizeof(double));
for (i = 0; i < data->n_sel; i++)
list[i] = data->sel[i]->attr;
mean = sample_mean(list, data->n_sel);
std = sample_std(list, mean, data->n_sel);
for (i = 0; i < data->n_sel; i++) {
d = (data->sel[i]->attr - mean)/std;
skewness += pow(d, 3.0);
kurtosis += pow(d, 4.0);
}
est[0] = skewness/data->n_sel;
est[1] = kurtosis/data->n_sel;
}
static double inverse_dist(DATA *data, DPOINT *where, double idPow) {
int i, j;
double sumweights, sumvals, value, weight, dist2;
static DATA *blockd = NULL;
if (data->n_sel <= 0)
ErrMsg(ER_IMPOSVAL, "zero neighbourhood in inverse_dist()");
if (data->n_sel == 1)
return data->sel[0]->attr;
blockd = block_discr(blockd, get_block_p(), where);
for (i = 0, value = 0.0; i < blockd->n_list; i++) {
for (j = 0, sumweights = sumvals = 0.0; j < data->n_sel; j++) {
dist2 = data->pp_norm2(data->sel[j], blockd->list[i]);
if (dist2 == 0.0) {
sumvals = data->sel[j]->attr;
sumweights = 1.0;
j = data->n_sel; /* break j loop */
} else {
if (idPow == 2.0)
weight = 1.0 / dist2;
else
weight = pow(dist2, -0.5 * idPow);
sumweights += weight;
sumvals += weight * data->sel[j]->attr;
}
}
value += blockd->list[i]->u.weight * sumvals/sumweights;
}
return value;
}
static void calc_poly(DPOINT *where, double *est) {
PLOT_POINT pt;
int j, k;
int n = get_n_edges();
POLYGON **edges = get_edges();
int *n_edges_polys = get_n_edges_polys();
pt.x = where->x;
pt.y = where->y;
for (j = 0; j < n; j++) { /* over edge files */
for (k=0; k<n_edges_polys[j]; k++) { /* over edges *//* for (j = 0; j < data->n_pol; j++) {*/
if (point_in_polygon(pt, &(edges[j][k]))) {
est[0]=j;
est[1]=k;
return;
}
}
}
}
void save_variogram_parameters(VARIOGRAM *v) {
int i;
if (vgm_pars == NULL)
vgm_pars = (double *) emalloc(3 * v->n_models * sizeof(double));
for (i = 0; i < v->n_models; i++) {
vgm_pars[3 * i + 0] = v->part[i].sill;
vgm_pars[3 * i + 1] = v->part[i].range[0];
vgm_pars[3 * i + 2] = v->part[i].range[1];
}
}
void reset_variogram_parameters(VARIOGRAM *v) {
int i;
for (i = 0; i < v->n_models; i++) {
v->part[i].sill = vgm_pars[3 * i + 0];
v->part[i].range[0] = vgm_pars[3 * i + 1];
v->part[i].range[1] = vgm_pars[3 * i + 2];
}
v->fit_is_singular = 0;
}